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1.  INTRODUCTION 


One  of  the  primary  tasks  in  the  present  effort  was  the  use  of 
the  S-CUBED  plastic-flow  code  to  calculate  for  typical  metal-matrix 
composites  tne  residual  microscopic  stresses  which  would  result  from 
given  fabrication  and/or  processing  cycles.  Such  calculations  have 
been  done  for  graphite/aluminum,  tungsten/aluminum  and  silicon- 
carbide/aluminum  uni  directionally  rein  jrced  fibrous  composites.  In 
particular,  a  rather  extensive  series  of  parametric  calculations  was 
performed  for  the  graphite/aluminum  case  to  determine  the  changes  in 
the  aluminum  matrix  residual  stresses  which  would  result  from  vari¬ 
ations  and/or  uncertainties  in  the  thermoelastic  properties  of  the 
graphite  fioers. 

Tne  S-CUBED  plastic-flow  code,  which  uses  the  concentric- 
cylinder  approximation  for  unidirectional  composites,  is  described 
in  detail  in  Ref.  1.  Briefly,  the  equations  of  elastic-plastic  flow 
are  integrated  numerically  as  the  composite  is  cooled  from  an 
initial  stress-free  state  and  then  reheated.  The  fibers  are  assumed 
to  remain  perfectly  elastic  throughout  the  cooling  and  heating 
cycle,  and  tney  may  be  treated  as  being  anisotropic- transversely 
isotropic.  The  matrix  material  is  assumed  to  be  isotropic  and  is 
allowed  to  undergo  plastic  flow  when  its  stress  state  reaches  the 
yield  surface.  The  magnitude  of  the  yield  stress  may  depend  both  on 
the  temperature  and  the  degree  of  plastic  flow. 

Results  of  the  parametric  study  of  a  graphite/aluminum-2024 
composite  are  given  in  Section  2,  where  stress-temperature  histories 
are  given  for  the  axial  stress,  hoop  stress,  and  fiber-matrix  inter¬ 
facial  radial  stress.  The  temperature  cycle  used  for  these  calcu¬ 
lations  is  a  cooling  from  a  consolidation  temperature  of  930°F  to  a 
minimum  temperature  of  -240°F  and  then  reheating  until  the  aluminum 
matrix  again  reaches  the  yield  surface.  Results  are  given  both  for 
the  TO  temper  (slow  cooling)  and  the  T4  temper  (rapid  quench) 
thermoelastic  properties  and  yield  strength  models  used  for  6061  and 
2024  aluminums  are  detailed  in  Appendix  A. 
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Results  of  two  calculations  on  a  tungsten/aluminum  composite 
are  given  in  Section  3,  and  similar  calculations  on  a  silicon 
carbide/aluminum  composite  are  described  in  Section  4. 

In  Section  5,  a  computational  method  is  presented  for  deter¬ 
mining  the  elastic  moduli  and  thermal  expansion  characteri sties  of 
cross-plied  laminates  from  the  properties  of  the  individual  uni- 
directionally  reinforced  composite  plies.  Calculated  results  as  a 
function  of  lay-up  angle  are  given  for  a  laminate  made  up  of 
graphite/aluminum  layers.  A  computer  code  has  been  constructed  for 
treating  a  laminate  with  plies  in  up  to  three  different  directions. 
The  properties  of  the  individual  plies  may  be  calculated  directly 
from  the  fiber  and  matrix  properties  if  desired.  The  code  also 
provides  values  for  tne  intra-ply  thermal  stress  derivatives  about 
an  initial  stress-free  state.  A  FORTRAN  listing  of  the  code  is 
given  in  Appendix  B. 

A  consideration  of  numerical  results  from  a  finite-element 
code  calculation  as  given  in  a  recent  report  by  Hashin  and  Humpnreys 
(Ref.  Z)  has  led  to  a  revision  of  the  PRUFC  code.  This  code  was 
recently  developed  at  S-CUBED  for  calculating  the  properties  of 
unidirectional  fibrous  composites  from  those  of  the  constituents 
(Ref.  3).  In  Section  6,  various  ways  of  obtaining  the  transverse 

properties  of  a  composite  within  the  context  of  the  concentric- 

cylinder  approximation  used  in  the  PRUFC  code  are  examined  and 
compared  with  the  results  of  a  more  nearly  exact  (but  much  more 
expensive  to  run)  finite-element  code  computation.  A  comparison  of 
residual  stresses  calculated  with  the  S-CUBED  plastic-flow  code 
(concentric-cylinder  approximation)  with  those  obtained  from  the 
finite-element  method  is  also  given  in  Section  6.  The  agreement  is 
quite  good. 

In  Section  7,  a  two-dimensional  Fourier-series  solution  is 

presented  for  the  stresses  and  thermal  deformations  in  an  elastic 
plate  subjected  to  an  arbitrary  transverse  temperature  di stribution. 

A  summary  and  some  general  conclusions  are  offered  in 

Section  8. 
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2.  GRAPHITE/ALUMINUM  PARAMETER  STUDY 


A  series  of  calculations  was  performed  to  determine  the 
sensitivity  of  the  matrix  residual  stresses  to  variations  in  various 
fiber  properties.  Eignteen  different  sets  of  graphite  fiber  pro¬ 
perties  were  used.  For  each  set  of  fiber  properties,  results  were 
ootained  for  a  slow  cooling  (TO  temper)  and  a  rapid  quench  (T4 
temper)  of  the  2024  aluminum  matrix  from  an  initial  consolidation 
temperature  of  930°F  to  a  minimum  temperature  of  -24Q°F.  The  fiber 
properties  are  taoulated  in  Table  2.1,  and  the  thenmoelastic  pro¬ 
perty  es  and  yield-stress  models  used  for  the  2024  aluminum  are 
sumnarized  in  Appendix  A.  The  room  temperature  thermoelastic 
properties  of  the  composite  as  calculated  witn  the  latest  version  of 
the  PRUFC  code  (Section  6)  for  each  set  of  fiber  properties  are 
given  in  Table  2.2,  the  values  being  appropriate  to  the  case  where 
the  aluminum  matrix  is  not  on  its  yield  surface. 

2.1  CALCULATIONAL  RESULTS 

Plots  of  tne  volume-averaged  (matrix)  axial  stress,  volume- 
averaged  hoop  stress,  and  the  radial  stress  at  the  fiber-matrix 
interface  as  a  function  of  temperature  are  given  as  Figures  2.1 
tnrough  2.42.  The  temperature  cycle  is  a  reheating  from  -240°F 
after  a  slow  (TO)  or  fast  (T4)  cooldown  from  an  initial  consoli¬ 
dation  temperature  of  930°F.  For  this  high  initial  temperature,  the 
matrix  is  completely  yielded  after  only  a  few  degrees  (<30°F)  of 
cooling,  and  kne  changes  in  the  stress  state  upon  further  cooling 
are  due  primarily  to  the  variation  of  yield  strength  with  temper¬ 
ature  and  plastic  flow.  Upon  reheating  from  -240°F,  however,  the 
matrix  immediately  drops  below  tne  yield  surface;  and  the  variation 
in  stress  with  temperature  is  then  governed  by  the  (reversible) 
equations  of  thermoelasticity.  As  the  reheating  continues,  the 
matrix  will  eventually  again  reach  its  yield  surface,  as  indicated 
oy  tne  open  circles  on  the  reheating  curves.  The  plots  are 
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TABLE  2.1 

FIBER  PROPERTIES  USED  IN  THE  PARAMETRIC  RESIDUAL-STRESS  STUDY 
OF  A  GRAPHITE/ALUMINUM  2024  UNIDIRECTIONALLY  REINFORCED  COMPOSITE 


ROOM- TEMPERATURE  G/AL  2024  COMPOSITE  PROPERTIES  AS  CALCULATED 
WITH  THE  PRUFC  CODE  (ELASTIC  MATRIX)  FOR  EACH  SET  OF 
FIBER  PROPERTIES  AS  GIVEN  IN  TABLE  2.1 


terminated  at  this  point  due  to  the  limitations  of  the  present 
vt  .ion  of  tne  S-CUBEO  plastic-flow  code.  Physically,  if  the 
reheating  were  continued  beyond  the  yield  point,  the  slope  of  the 
stress-temperature  curves  would  exhiDit  a  discontinuous  sharp 
decrease  in  absolute  magnitude  with  a  subsequent  return  to  zero 
stress  at  the  melt  point.  As  long  as  tne  temperature  at  tne  upper 
yield  point  is  not  exceeded,  subsequent  cooling  will  cause  the 
stresses  to  retrace  the  original  heating  path  from  -240°F.  (The 
possible  effect  of  creep  on  the  matrix  residual  stresses  is  not 
treated  nere.) 

Tne  effect  on  tne  matrix  stress  of  a  variation  of  tne  axial 
coefficient  of  thermal  expansion  for  the  fibers  is  shown  in  Figures 
2.1  through  2.6.  Meitner  the  axial,  hoop,  nor  interfacial  radial 
stress  snows  a  marked  variation  for  the  range  a  =  0.5  x  1Q~^ 

I  C  1  ^ 

°F“  to  a  *  -  1.5  x  10  F~  for  the  fiber  expansion  coef¬ 
ficient.  The  reversiule,  elastic  portion  of  the  axial  stress  shows 
the  greatest  dependency. 

Tne  effect  of  a  variation  in  the  fibers  transverse  coefficient 
of  thermal  expansion  on  the  matrix  stresses  is  depicted  in  Figures 
2.7  tnrough  2.12.  The  axial  stress  does  not  exhibit  a  strong  depen¬ 
dency,  but  the  hoop  and  interfacial  radial  stresses  show  a  marked 
variation.  Perhaps  it  is  not  surprising  that  as  the  fiber's  trans¬ 
verse  expansion  coefficient  approaches  that  of  the  matrix,  the  hoop 
and  radial  residual  stresses  become  small.  Here  it  should  be 
pointed  out  that  a  negative  (compressive)  interfacial  radial  stress 
would  increase  the  effective  fiber-matrix  bond  strengtn;  hence  a 
moderately  large  negative  value  might  be  desirable. 

A  variation  of  tne  fiber's  Poisson's  ratios,  as  shown  in 
Figures  2.13  tnrough  2.18,  has  only  a  minute  effect  on  all  three  of 
the  calculated  stress-temperature  histories. 

An  increase  in  the  fiber's  axial  elastic  modulus  from  50  Msi 
to  100  Msi,  Figures  2.19  through  2.24,  does  produce  a  more  rapid 
decrease  in  the  matrix  axial  stress  upon  reheating  from  -240°F  but 
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has  a  completely  negligiDle  effect  on  the  hoop  and  interfacial 
radial  stress-temperature  histories. 

A  parametric  variation  in  the  transverse  elastic  modulus  of 

the  fibers  leads  to  the  calculated  results  for  the  matrix  stresses 
plotted  in  figures  2.25  througn  2.30.  Doubling  the  fiber  modulus 
from  1  Msi  to  2  Msi  has  a  negligible  effect  on  the  axial  matrix 
stresses,  and  the  effect  on  the  calculated  hoop  and  interfacial 
radial  stresses  is  also  small. 

For  tne  calculations  discussed  above,  the  fiber  volume  frac¬ 
tion  in  all  cases  is  n^  =  0.45.  In  Figures  2.31  through  2.36, 
results  are  plotted  for  a  fiber  volume  fraction  of  n^  =  0.37  for 
two  values  of  the  fiber's  transverse  thermal  expansion  coefficient. 
As  was  tne  case  with  Calculations  2,  4  and  6,  the  larger  value  of 
tne  transverse  expansion  coefficient  produces  a  significant  decrease 
in  tne  magnitudes  of  the  hoop  stresses  and  interfacial  radial 

stresses. 

Results  for  two  different  fiber  volume  fractions  of  nf  = 
0.30  and  n^  =  0.60  are  plotted  in  Figures  2.37  through  2.42.  The 
residual  stresses  are  somewhat  higher  in  magnitude  for  the  larger 
volume  fraction. 

2.2  SUMMARY  AMD  C0MCLUSI0MS 

In  general,  for  all  of  the  fiber  parameter  variations  dis¬ 
cussed  aoove,  the  calculations  indicate  that  the  axial  matrix  stress 
at  roan  temperature  after  a  rapid  quencn  from  930°F  to  -240°F  (T4 

temper)  will  lie  in  the  range  from  about  30  ksi  to  40  ksi  (tensile) 

and  tnat  tne  matrix  will  remain  elastic  up  to  about  600°F  upon 
reheating.  The  residual  hoop  stress  at  room  temperature  in  all 
cases  is  tensile  and  ranges  from  about  10  ksi  to  20  ksi  in  magnitude 

if  the  fiber's  transverse  coefficient  of  thermal  expansion  is  taken 

-6°  -1 

to  be  less  than  or  equal  to  that  of  the  matrix  ( ~  1 2 . 4  x  10  F 

for  2024  aluminum  at  room  temperature) .  Further,  the  interfacial 
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radial  stress  Detween  the  fiber  and  matrix  is  compressi ve  as  long  as 
the  fioer's  transverse  expansion  coefficient  does  not  exceed  tnat  of 
the  matrix,  the  magnitude  ranging  from  about  2  ksi  to  10  ksi  over 
tne  range  of  material  parameters  used  in  these  calculations.  This 
would  appear  to  be  a  significant  result,  since  a  compressive  inter- 
facial  stress  is  desirable  as  an  aid  in  maintaining  the  integrity  of 
tne  fiber-matrix  bond.  This  Dond  is  critical,  since  it  is  the 
source  of  tne  (localized)  snear  forces  which  a^e  ultimately  respon¬ 
sible  for  equilibrating  axial  strain  between  the  fiber  and  matrix 
components. 

For  tne  slow  cool  to  -240°F  (To  temper)  the  axial  residual 

stresses  at  room  temperature  for  all  cases  considered  here  are  very 
close  to  zero,  a  situation  which  would  appear  to  be  desiraole  if  the 
composite  is  to  be  subjected  to  subsequent  mechanical  loading. 

Because  of  tne  reduced  yield  strength  of  the  TO  temper  with  respect 
to  that  of  the  T4,  however,  tne  matrix  will  again  reach  its  yield 
point  after  reneating  to  between  200°F  and  3J0°F.  Even  tnough  the 
axial  stress  at  room  cemperature  is  much  lower  for  the  TO  temper 

than  tne  T4,  the  hoop  stress  will  still  De  dDout  5  ksi  to  10  Ksi  in 

tension  for  the  range  of  fiber  parameters  considered.  As  long  as 
tne  fiber's  transverse  expansion  coefficient  is  less  than  tnat  of 
tne  matrix,  the  interfacial  radial  stresses  would  still  be 
compressive  according  to  the  calculations,  the  magnitudes  being  in 
the  range  from  2  ksi  to  about  8  ksi. 

A  room  temperature  axia1  residual  stress  of  30  ksi  to  40  ksi 
in  the  matrix  as  calculated  for  the  T4  temper  would  be  undesirable 
if  tne  composite  were  to  be  subjected  to  large  tensile  loads.  This 
is  due  to  the  fact  that  the  matrix  would  reach  the  yield  surface, 
witn  a  consequent  degradation  of  the  axial  modulus,  at  a  lower 
applied  load  tnan  it  would  if  the  matrix  stress  were  initially  near 
zero.  For  compressive  loads  about  the  ambient  state,  however,  an 
initial  tensile  residual  stress  might  be  useful.  In  any  case,  it  is 
not  immediately  apparent  how  the  room  temperature  residual  stress 
for  the  T4  temper  might  be  reduced.  Continued  heating  beyond  the 
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yield  point  of  ~600°F  would  result  in  a  stress  temperature  path  that 
must  ultimately  pass  through  zero  stress  at  the  melt  point.  A 
recooling  from  any  point  on  this  path  segment  would  result  in  a 
(reversible)  elastic  curve  for  the  axial  stress  that  would  lie  to 
the  right  of  the  initial  heating  curve  from  -24G°F,  and  consequently 
the  room  temperature  axial  residual  stress  would  be  higher  than  it 
was  before.  This  might  oe  alleviated  however,  if  creep  of  the 
matrix  material  was  taken  into  account. 

Some  reduction  of  the  axial  residual  stress  would  be  achieved 
oy  couling  below  -240°F  initially.  No  explicit  calculation  was 
done,  but  it  appears  from  an  extrapolation  of  the  curves  presented 
nere  that  cooling  to  ~-S20°F  (liquid  nitrogen)  would  reduce  tne 
axial  residual  stress  by  about  5  ksi  or  10  ksi . 
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AXIAL  STRESS  (ksi) 
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Temperature  (°F) 

Figure  2.2.  Volume-averaged  axial  stress  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite  (slow 
cool ing) . 


11 


</l 


c/1 
c r> 

LlJ 

cc. 
I — 

C/1 

Q- 

O 

O 


0  500  1000 


Temperature  (°F) 


Figure  2.3.  Volume-averaged  hoop  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid  cooling). 
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HOOP  STRESS  (ksi ) 


Temperature  (°F) 

Figure  2.4.  Volume-averaged  hoop  stress  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite  (slow 
cool  i  ng) . 
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Temperature  (°F) 

Figure  2.5.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphite/aluminum  composite  (rapid 
cool ing) . 


14 


Temperature  (°F) 

Figure  2.6.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphi te/aluminum  composite  (slow 
cool i ng) . 
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Figure  2.7.  Volume-averaged  axial  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid 
cool  ing) . 
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AXIAL  STRESS  (ksi) 


Figure  2.8. 


Volume-averaged  ax^al  sH;resi  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite  (slow 
cool i ng) . 
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Figure  2.10.  Volume-averaged  hoop  stress  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite  (slow 
cool ing) . 
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INTERFACIAL  RADIAL  STRESS  (ksi) 
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Figure  2.12.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphite/aluminum  composite  (slow 
cool i ng) . 
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Figure  2.15.  Volume-averaged  hoop  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid 
cool i ng) . 
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Figure  2.16.  Volume-averaged  hoop  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (slow 
cool i ng) . 


25 


Temperature  (°F) 

Figure  2.18.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  qraphi te/aluminum  composite  (slow 
cool ing) . 
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Figure  2.19.  Volume-averaged  axial  stress  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite 
(rapid  cool ing) . 
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Figure  2.20.  Volume-averaged  axial  stress  versus  temperature 
for  the  matrix  of  a  graphite/aluminum  composite 
(slow  cool i ng) . 
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Figure  2.21.  Volume-averaged  hoop  stress  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite  (rapid 
cooling) . 
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Figure  2.22.  Volume-averaged  hoop  stress  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite 
(slow  cool ing) . 
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Figure  2.23.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  *  graphi te/aluminum  composite 
(rapid  cool ing) . 
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Figure  2.25.  Volume-averaged  axial  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid  cooling). 
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Figure  2.27.  Volume-averaged  hoop  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid 
cool ing) . 
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Figure  2.28. 
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Figure  2.29.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphi te/al umi num  composite  (rapid 
cool i ng) . 


INTERFACIAL  RADIAl  STRESS  (ksi) 


AXIAL  STRESS  (ksi) 


Figure  2.31.  Vol ume-averaged  axial  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid  cooling). 


Figure  2.32.  Volume-averaged  axial  stress  versus  temperature  for 
the  matrix  of  a  graphite/aluminum  composite  (slow 
cool ing) . 


41 


Temperature  (°F) 

Figure  2.33.  Volume-averaged  hoop  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid  cooling). 
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Figure  2.34.  Volume-averaged  hoop  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (slow  cooling). 
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Figure  2.35.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphite/aluminum  composite  (rapid 
cooling) . 


Figure  2.36.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphite/aluminum  composite  (slow 
cooling) . 
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AXIAL  STRESS  (ksi) 


Figure  2.37.  Volume-averaged  axial  stress  versus  temperature  for  the 
matrix  of  a  graphite/aluminum  composite  (rapid  cooling). 
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Figure  2.38.  Volume-averaged  axial  stress  versus  temperature  for  the 
matrix  of  a  graph i te/al umi num  composite  (slow  cooling). 
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Figure  2.41.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphite/aluminum  composite  (rapid 
cooling) . 
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Fiyure  2.42.  Radial  stress  versus  temperature  at  the  fiber-matrix 
interface  of  a  graphite/aluminum  composite  (slow 
cool ing) . 
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3.  TUNGSTEN/ ALUMINUM  RESIDUAL  STRESSES 


Calculated  axial  and  hoop  stresses  (volume  averaged)  in  the 
aluminum  matrix  are  given  in  Figures  3.1  and  3.2  for  a  composite 
whicn  is  45  percent  by  volume  of  tungsten  fibers  and  in  Figures  3.3 
and  3.4  for  a  55  percent  fiber  loading.  The  temperature  cycle  is  a 
cooling  from  a  consolidation  temperature  of  930°F  to  a  minimum 
temperature  of  -240 °F  followed  by  a  reheating.  The  open  circles  on 
the  reneating  curves  mark  the  temperature  at  which  the  stress  state 
in  tne  matrix  again  reaches  the  yield  surface.  As  was  the  case  for 
tne  graphite/aluminum  calculations  in  Section  2,  the  matrix  reaches 
the  yield  surface  after  only  a  few  degrees  of  cooling  (<  50°F)  from 
the  initial  value  of  930°F.  The  suDsequent  increase  in  stress  is 
due  to  the  temperature  and  plastic-flow  dependence  of  the  yield 
strength. 

The  material  properties  and  yield-strength  models  for  the  2024 
aluminum  are  given  in  Appendix  A.  The  tungsten  fibers  are  assumed 
to  be  elastic  over  tne  entire  temperature  range,  the  following 
values  being  used  for  the  tungsten's  elastic  modulus,  Poisson's 
ratio  and  linear  coefficient  of  thermal  expansion: 

E  -  59  Msi 

v  =  0.2b 

a  =  2.5  x  10"6  in/in/°F  . 

Tne  properties  of  the  composite  as  calculated  with  the  PRUFC  code 
are  given  in  Table  3.1  below. 
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TABLE  3.1  CALCULATED  THERMOELASTIC  PROPERTIES 
OF  A  TUrfGSTEM/AL2024  COMPOSITE 


nf 

0.45 

0.55 

£x(Msi ) 

17.74 

20.69 

£2(Msi ) 

32.33 

37.18 

vxy 

V  _ 

0.434 

0.428 

0.305 

0.300 

.WV1) 

8.48 

7.34 

az(10'6*F'1) 

4.34 

3.82 

axy(Msi) 

Gx2(Msi) 

6.19 

7.24 

7.63 

8.98 

For  the  TO  temper  {slow  cooling  from  930°F)  the  residual  axial 
and  hoop  stresses  are  close  to  zero  after  reheating  to  room  tempera¬ 
ture  from  -240°F.  Because  of  the  low  yield  strength,  however,  the 
matrix  will  again  reach  the  yield  surface  at  about  150°F.  For  the 
T4  temper  (rapid  quench),  the  higher  yield  strength  allows  the 
matrix  to  remain  elastic  up  to  about  500°F.  The  high  value  of  the 
axial  residual  stress  at  room  temperature,  however,  would  be 
undesirable  if  the  composite  is  to  be  subjected  to  large  tensile 
1  oads. 
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Figure  3.2.  Volume-averaged  axial  and  hoop  stresses  versus 

temperature  for  the  matrix  of  a  tungsten/aluminum 
composite  (slow  cooling). 
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0  500  1000 

Temperature  (°F) 

Figure  3.3.  Volume-averaged  axial  and  hoop  stresses  versus 

temperature  for  the  matrix  of  a  tungsten/aluminum 
composite  (ranid  cooling). 
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4.  SILICON  CARBIDE/ALUMINUM  RESIDUAL  STRESSES 


Several  calculations  have  been  done  to  determine  the  magni¬ 
tudes  of  the  residual  stresses  to  be  expected  in  a  SiC/Al  composite 
as  a  result  of  cooling  from  the  consolidation  to  ambient  tempera¬ 
ture.  The  effect  of  these  stresses  on  the  elastic  properties  of 
unidirectional  panels  and  cross-plied  laminates  is  also  considered. 
The  results  of  this  section  have  been  submitted  previously  as  an 

Interim  Progress  Report  to  the  Naval  Research  Laboratory  (Ref. 4). 

4.1  FIBER  PROPERTIES 

For  bulk  SiC,  whicn  we  assume  to  be  isotropic,  the  elastic 
(Young's)  modulus  and  shear  modulus  were  taken  as  (Ref.  5): 

£  =  68  Msi 
S  =  27  Msi  . 

Tne  coefficient  of  tnermal  expansion  was  assumed  to  be 
a  =  2.42  x  10"6  in/i n/°F“1 

which,  from  Ref  5,  is  appropriate  for  the  temperature  range  con¬ 

sidered  here. 

The  SiC  fibers,  however,  are  not  completely  bulk  material,  but 

apparently  contain  a  core  of  essentially  zero  strength  material 

-3  . 

(carbon)  whose  diameter  is  about  1.3  x  10  in.  as  compared  with 

_  3 

an  outer  fiber  diameter  of  about  5.6  x  10  in.  The  volume 
fraction  of  the  caroon  core  is  thus  0.054.  The  overall  properties 
of  a  fiber  were  determined  with  the  PRUFC  Code  (Ref.  3),  in  which 
tne  sheatn  of  the  concentric  cylinder  model  was  assigned  properties 
appropriate  to  the  bulk  SiC,  and  the  core  was  assigned  elastic 

moduli  a  factor  of  10  smaller  (hollow  cylinder  approximation) . 
The  resulting  calculated  properties  of  the  fiber,  which  now  becomes 
ani sotropic-transversely  isotropic,  are: 
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ea 

=  64.33  Msi 

£t 

=  60.56  Msi 

VTT 

=  0.2400 

VTA 

=  0.2600 

gat 

=  24.25  Msi 

aA  =  aT  =  2.42  x  lO'^F'1  , 

where  the  Poisson's  ratio  v-j-j  refers  to  the  transverse  contraction 
for  an  applied  transverse  stress,  and  vjA  refers  to  the  transverse 
contraction  for  an  applied  axial  stress.  As  expected,  because  of 

tne  small  volume  fraction  of  the  core,  the  net  fiber  properties 

differ  only  slightly  from  those  of  the  bulk  SiC. 

Tne  material  properties  and  yield-strength  models  for  the 
aluminum  6061  matrix  are  given  in  Appendix  A. 

4.2  RESIDUAL  STRESS  RESULTS 

Average  axial  stresses  in  the  aluminum  matrix  of  a  SiC/Al 
composite  are  shown  in  Figure  4.1.  Results  are  given  for  fiber 
volume  fractions  of  0.35  and  0.50  and  for  both  the  TO  and  T4 

aluminum  tempers.  The  assumed  thermal  cycle  process  is  a  cooling 

from  an  initially  stress-free  state  at  700°F  to  room  temperature  and 
then  reheating. 

As  mentioned  previously,  the  code  uses  a  concentric-cylinder 
approximation  to  model  the  composite.  For  initial  cooling  from  the 
stress-free  state,  both  the  inner  cylinder  (fiber)  and  the  outer 
(matrix)  are  elastic.  As  cooling  continues,  tne  matrix  material 
will  in  general  reach  its  yield  point,  after  which  the  Reuss 
equations  for  plastic  flow  (see,  for  example.  Ref.  6)  are  integrated 
numerically  to  determine  the  distribution  of  stress  in  the  outer 
(matrix)  cylinder.  In  the  code's  present  form  the  fiber  is  treated 
ii.->oariy  elastic  over  the  entire  temperature  range. 
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Figure  4.1.  Calculated  average  axial  stresses  in  the  matrix  of  a  SiC/Al  6061  composite  as  a 
function  of  temperature  for  different  fiber  volume  fractions  and  matrix  tempers. 
The  assumed  thermal  cycle  is  a  cooling  to  70°F  from  a  consolidation  temperature 
of  700°F,  and  then  a  reheating  from  70°F  to  indicate  the  rate  of  decrease  of  the 
residual  stress. 


For  the  cases  depicted  in  Figure  4.1,  about  20  degrees  or  less 
of  cooling  from  the  initial  temperature  of  700°F  is  sufficient  to 
bring  the  aluminum  matrix  to  its  yield  surface.  The  additional 
increase  in  axial  stress  for  further  cooling  results  from  two 
causes:  the  increase  of  tne  matrix  yield  stress  itself  as  a  func¬ 
tion  of  decreasing  temperature,  and  a  further  increase  due  to 
plastic-flow  work  hardening.  At  the  lowest  temperature  { 70 °F )  the 
axial  stress  is  a  maximum,  and  its  magnitude  depends  strongly  on  the 
effective  temper  of  the  aluminum.  The  calculation  suggests  that  the 
maximum  residual  stress  due  to  rapid  quencning  (T4  temper)  would  be 
over  twice  that  of  the  slow-cooled  or  annealed  state.  If  the  com¬ 
posite  is  reheated  from  room  temperature,  the  matrix  material  will 
fall  off  the  yield  surface,  and  as  shown  in  Figure  4.1,  the  rate  of 
decrease  of  stress  in  the  elastic  regime  will  be  quite  rapid.  For 
the  annealed  case,  a  temperature  increase  of  about  100°F  would  be 
sufficient  to  reduce  the  residual  stresses  to  approximately  zero. 

It  is  anticipated  that  A357  aluminum  casting  alloy  will  also 
be  used  in  the  fabrication  of  these  composites.  No  residual  stress 
calculations  for  this  material  have  been  done  as  yet  jecause  we  have 
been  unable  to  locate  yield  strength  information  for  the  TO  and  T4 
tempers.  No  large  differences  witn  respect  to  the  6061  matrix  are 
expected  however,  since  the  yield  strengtn  of  the  A357  alloy  in  the 
T6  temper  is  comparable  to  that  of  6061-T6,  and  other  mechanical 
properties  also  are  not  signi ficantly  different. 

4.3  ELASTIC  PROPERTIES  OF  A  UNIAXIALLY  REINFORCED  LAMINATE 

Since  tne  coefficient  of  thermal  expansion  for  the  aluminum 
matrix  is  greater  than  that  of  the  Si C  fibers,  the  residual  stresses 
in  tne  matrix  after  cool-down  will  be  tensile  in  the  axial  direc¬ 
tion.  (The  calculated  hoop  stresses  in  the  concentric  cylinder 
model  are  also  tensile,  but  the  radial  stress  is  compressive,  so 
that  there  should  be  little  tendency  for  deoonding  to  occur  between 
tne  matrix  and  fibers.)  If  the  composite  is  not  subjected  to  an 


overall  tensile  axial  stress,  the  matrix  will  remain  on  the  yield 
surface  and  the  effective  axial  elastic  modulus  will  be  lower  than 
it  would  be  in  the  absence  of  residual  stresses. 

The  plastic-flow  code  as  described  in  Ref.  1  has  been 
augmented  so  that  it  is  now  possible  to  calculate  the  stresses  as  a 
function  of  axial  strain  following  the  cool-down  procedure.  The 
calculated  overall  axial  stress  as  a  function  of  axial  strain  at 
70  F  is  plotted  in  Figure  4.2  for  tne  6061-T0  yield  model  after 
cooling  from  700°F.  Also  shown  is  the  straight-line  stress-strain 
curve  for  an  elastic  matrix  as  calculated  with  the  PRUFC  code,  and 
an  experimental  stress-strain  curve  obtained  from  AVCO  (Ref.  7),  the 
volume  fraction  of  fibers  in  this  case  being  0.45.  The  experimental 
curve  changes  slope  at  about  90  ksi,  and  the  slope  of  the  higher- 
stress  segment  is  in  good  agreement  with  that  of  the  curve  corre¬ 
sponding  to  the  yielded  matrix.  The  lower-stress  portion  of  the 
experimental  curve  however,  is  in  agreement  with  the  stress-strain 
curve  as  calculated  using  an  elastic  matrix. 

The  overall  axial  stress  for  the  composite  is  related  to  that 
in  the  fibers  and  matrix  by 

°z  =  Vzf  +  vi  '  nf)azm  * 

For  this  material  one  has,  to  good  approximation, 

.  Ef 
^  •» 

for  tne  fiber  stress  in  terms  of  that  in  tne  matrix,  where  E^  and 
Em  are  the  respective  fiber  and  matrix  axial  elastic  moduli.  For 
an  elastic  matrix  of  modulus  Em  =  10  Msi,  and  an  overall  composite 
axial  stress  of  90  ksi,  the  matrix  would  experience  a  stress  of  26 
ksi,  a  value  close  to  the  nominal  yield  stress  of  6061  aluminum  in 
the  T4  temper.  On  tne  other  hand,  if  the  Dreak  in  the  experimental 
curve  at  90  ksi  is  to  be  interpreted  as  the  transition  point  from  an 
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Comparison  of  calculated  and  experimental  axial  stress-strain  curves  for  a 
0.45  fiber  volume-fraction  SiC/Al  6061  composite.  The  curve  for  the  yielded 
matrix  corresponds  to  the  TO  temper  model  of  Figure  4.1. 


elastic  to  a  plastic  matrix,  then  the  initial  state  of  the  matrix 
would  have  to  be  essentially  stress-free,  in  marked  disagreement 
with  the  calculated  axial  residual  stress  as  shown  in  Figure  4.1. 

The  slope  of  the  calculated  stress-strain  curve  corresponding 
to  a  yielded  matrix  is  29.3  Msi,  which  is  close  to  the  value 

nf£f  =  (0.45) (64.33)  =  28.9  Msi  , 

which  one  would  expect  for  the  composite  elastic  modulus  if  the 
matrix  modulus  were  zero.  This  result  is  to  be  expected,  since  at 
tne  yield  surface  the  effective  modulus  of  the  matrix  is  approxi¬ 
mately  equal  to  the  local  slope  of  the  yield-strength  versus  strain 
curve,  wnich  typically  ranges  from  several  hundred  ksi  to  about  1 
Msi  (see  Figures  A. 3  and  A. 4  of  Appendix  A).  Since  this  slope  for 
either  tne  TO  or  T4  temper  is  small  as  compared  with  the  64  Msi 
modulus  of  the  fiber,  the  choice  of  tempers  has  a  negligible  effect 
on  tne  calculated  value  of  the  elective  axial  modulus  for  the 
composite  if  the  matrix  is  indeed  on  the  yield  surface. 

The  concentric-cylinder  approximation  used  in  the  S-CUBED 
plastic-flow  code  does  not  lend  itself  to  the  calculation  of  the 
transverse  stress-strain  relations  if  the  matrix  is  on  the  yield 
surface.  From  the  PRUFC  code,  however,  for  an  elastic  matrix  the 
calculated  transverse  modulus  for  the  composite  is 

Er  =  17.1  Msi 


for  a  fiber  volume  fraction  of  0.45.  The  mixture  result  for  a 
simple  series  combination  of  tne  fioer  and  matrix  materials. 


0.45  . 
60.56 


0.55 

~T0~" 


gives  the  result 


£f  =  16.0  Msi  , 


quite  close  to  the  PRUFC  result,  so  that  the  simple  series  combina¬ 
tion  might  be  expected  to  give  a  useful  approximation  to  the 
transverse  behavior,  even  if  the  matrix  has  yielded. 

Note  here  that  if  the  matrix  material  is  at  the  yield  surface 
due  to  trie  cool-down  residual  stresses,  then  the  dominant  residual 
stress  is  tension  in  the  axial  direction.  Consequently,  an  applied 
transverse  tension  will  pull  tne  matrix  off  the  yield  surface  and 
the  initial  transverse  behavior  will  be  elastic.  Under  an  applied 
transverse  compressive  stress,  the  matrix  would  stay  on  the  yield 
surface  and  flow  plastically,  resulting  in  a  lower  apparent  modulus 
for  the  composite. 

For  the  simple  series  combination,  the  overall  strain  is  given 
by 

e  =  nfef  +  (i  "  nf,em 


in  the  transverse  direction,  or 


E  =  nf^M1-  nf)em(a) 


The  above  relation  was  used  to  obtain  the  estimates  for  the  90 
degree  orientation  composites  as  plotted  in  Figure  4.2. 


4.4  LAMINATE  PROPERTIES 

The  code  constructed  at  S-CU8ED  (Appendix  B)  may  be  used  to 
compute  the  elastic  properties  of  a  laminate  with  plies  oriented  in 
up  to  three  different  arbitrary  directions.  As  input  one  may  use 
either  the  elastic  and  thermal  properties  of  each  type  of  ply,  or 
tne  volume  fractions  and  individual  properties  of  the  fiber  and 
matrix  components  for  each  type  of  ply.  In  the  latter  case  the  code 
computes  the  ply  properties  as  per  the  method  used  in  the  PRUFC 
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code.  With  the  fiber  and  matrix  properties  as  input,  the  code  was 
used  to  determine  the  elastic  behavior  of  a  0°/90o/90o/0°  laminate 
with  each  ply  consisting  of  a  0.45  fiber  volume-fraction  SiC/Al 
composite.  The  resulting  slope  of  tne  initial  elastic  segment  is 
indicated  in  Figure  4.3,  both  for  an  initial  stress-free  state  and 

for  tne  case  where  tne  matrix  material  is  at  the  yield  surface  due 

to  cool-down  residual  stresses.  Fcr  the  90°  plies,  the  initial 
benavior,  as  discussed  above,  will  be  elastic  even  tnough  the  matrix 
is  initially  on  the  yield  surface.  It  turns  out  that  for  these 

orientations,  tne  resulting  elastic  modulus  of  the  composite 
laminate  is  very  nearly  the  same  as  that  of  a  simple  combination  of 
parallel  springs.  Thus,  with  equal  volume  fractions  of  the  0°  and 
90°  plies,  the  laminate  modulus  correspondi ng  to  an  elastic  matrix  is 

(34.5  +  17.1)/2  =  25.8  Msi, 

and  for  the  yielded  matrix  it  drops  to 

(29.3  +  17.1J/2  =  23.2  Msi 

as  indicated  in  Figure  4.3.  Eventually,  the  90°  plies,  which  had 

initially  dropped  off  the  yield  surface,  will  once  again  reach  the 
yield  surface.  Using  an  elastic-perfectly  plastic  model  with  a 
yield  stress  of  21  k si  (T4  temper),  this  is  found  to  occur  at  an 
overall  stress  of  aoout  28  ksi ,  at  which  point  the  slope  of  tne 
stress-strai n  curve  will  drop  to  about 

29.3/2  =  14.6  ksi 

Tne  resulting  stress-strai n  relationship  agrees  well  with  a 
projected  estimate  from  AVCO,  as  shown  in  Figure  4.3. 


250. 
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Figure  4.3.  Calculated  stress-strain  relations  for  a  0°/90o/90o/0°  laminate  consisting  of 
0.45  fiber  volume-fraction  SiC/Al  plies.  Also  shown  is  a  projected  estimate 
from  AVCO  (Ref.  7). 


5.  THERMOELASTIC  PROPERTIES  OF  CROSS-PLY  LAMINATES 


Computational  machinery  is  presented  for  determining  the 
elastic  moduli  and  thermal  expansion  characteristics  of  cross-plied 
laminates  from  the  properties  of  the  individual  uni  directionally 
reinforced  layers.  Explicit  formulas  are  given  for  the  elastic 
constants  and  thermal  expansion  coefficients  for  a  laminate  with 
equal  volume  fractions  of  plies  oriented  in  two  directions.  The 
computer  code  (Appendix  B)  in  its  present  form  may  be  used  for  plies 
oriented  in  up  to  three  different  directions.  The  code  also  pro¬ 
vides  the  stress-temperature  derivatives  within  each  type  of  ply, 
appropriate  to  a  uniform  heating  from  an  initial  stress-free  state. 

5.1  TRANSFORMATION  OF  ELASTIC  CONSTANTS 

As  a  first  step,  it  is  necessary  to  obtain  expressions  for  the 
elastic  constants  of  a  fiber-reinforced  material  in  a  coordinate 
system  rotated  with  respect  to  that  defined  by  the  direction  of  the 
fibers.  The  stress-strain  relations  for  a  f iber-reinforced  material 
assume  the  simplest  form  in  the  fiber  coordinate  system  (FCS),  i.e., 
tnat  coordinate  system  in  which  one  of  the  axes  is  in  the  direction 
of  the  fibers.  If  we  take  the  z-axis  as  being  parallel  to  the 
fibers,  then  tne  stresses  are  given  in  terms  of  the  strains  by  the 
matrix  relation 
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"xx  =  Cllexx  *  ^12eyy  *  C13szz 
°yz  *  2C44eyz  ’  etc" 

where  the  strains  are  defined  by 


(5.2) 


(5.3) 


In  what  follows,  the  subscripts  1,  2,  3  will  be  used  interchangeably 
with  x,  y,  z,  respecti vely ,  to  denote  the  coordinate  axes.  In  the 
usual  case,  where  the  material  is  not  only  orthotropic  but  also 
transversely  isotropic, 

C66  =  I  *C11  "  C12^  *  ^5*4* 


For  tne  laminate  analysis,  we  require  the  elastic  constants  of 
the  ply  in  the  rotated  coordinate  system  x'y'z'  as  shown  in  Figure 
5.1,  obtained  by  a  positive  rotation  (e  >  0)  about  the  x,x'  axis, 
which  is  perpendicular  to  the  face  of  the  ply.  The  elastic  con¬ 
stants  C  are,  of  course,  components  of  the  fourtn  rank  tensor, 
pq 

v'ijmn*  which  reldtes  the  stresses  to  the  strains, 

a . •  =  C . .  e  .  (5.5) 

ij  ijmn  mn 


using  the  summation  convention  m,n  =  1,2,3  over  repeated  indices. 


In  general,  the  transformation  to  a  new  coordinate  system  is 
given  by 


C 


_  a  off  f  r 

pqkr  ~  pi *qj  km  rn  ijmn 


(5.6) 


where  the  2...  relate  the  components  of  a  vector  in  the  primed 

'  J 

coordinate  system  to  those  in  the  unprimed;  and  in  the  present  case 
(Figure  5.1),  these  quantities  are  given  by  the  matrix 
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—  pn  hi 

£ij  =  hi  hi 

0  0 

J  |_  31  32 

In  the  FCS,  the  nonvanishing  components  of  the  tensor  C. 
are  given  in  terms  of  the  elastic  constants  defined  in  Eq,  (5.1)  by 

C11  =  CU11  =  C2222 

C12  =  C 1 122  =  C2211 

C13  =  C1133  =  C2233  =  C3311  =  C3322 

r  „  (5.8) 

u33  =  ^3333 

C44  =  C2323  =  C2332  =  C32 23  =  C3232 

=  C1313  =  C1331  =  C31 13  =  C3131 

C66  =  C1212  =  C1221  =  C21 12  C2121  * 

Thus,  tne  calculation  of  a  typical  component  in  the  new  coordinate 
system  requires  the  consideration  of  21  terms.  For  example, 

C3333  =  £31£31£31£31C1111  +  £32£32£32£32C2222 

+  £31£31£32£32C1122  +  £32£32£31£31  C2211  (5.9) 

+  #o<  +  . . .  + 

*  9  9  9  9  r  +  9  9  o  9  r* 

32  31  31  32^2112  *32  31  32  31^2121  * 

Many  of  the  terms  are  zero,  however,  because  of  the  simple  form  of 
the  transformation  matrix,  Eq.  (5.7). 

The  result  is  that  the  stresses  and  strains  in  the  primed 

coordinate  system  of  Figure  5.1  are  related  by  (note  that  the  matrix 

is  not  symmetric) 
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where  the  new  elastic  constants  are  given  in  terms  of  the  old  by 


ch 

■Cll 

CM 

H 

o 

II 

Ch 

C12 

Ci4 

=  -  C 

C22 

=  cn< 

>2e  +  C13sin2( 


2  2 

+  4C44sin  ecos  e 

2  2  4  4 

0^3  =  (C31  +  C33)  sin  ecos  e  +  C13(sin  e  +  cos  e) 

2  2 

-  4C44sin  ecos  e 

3  3  3 

C24  =  ”  cns'inecos  9  +  C^3(  sinecos  e  -  sin  ecose) 

3  3  3 

+  C33sin  ecose  +  2C^(  sinecos  e  -  sin  eco^e) 

4  4  2  2 

C33  =  C33cos  e  +  C^sin  e  +  2C^3sin  ecos  e 

2  2 

+  4C44sin  9005  9 

3  3 

C34  =  "  cnsl'n  90059  +  C33sinecos  e 

3  3 

+  C  13(  si n  ecose  -  sinecos  e) 

3  3 

♦  2C^(sin  ecose  -  sinecos  e) 


(5.11) 
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u44  = 


C'  - 

33 


°56 

C6o 


2  2 

(Cn  +  C^-j)  sin  ecos  e 

+  C44(cos2e  -  sin2e)2 

2  2 
C^cos  e  +  C66sin  e 

(C44  '  C6o)  sin9C0S9 

C66cos2®  +  C44sin2® 


2  2 
2C13sin  ecos  e 


The  above  results  have  been  checked  against  tnose  given  by 

r  8 1 

HasninL  for  an  initial  coordinate  system  in  wnich  the  fibers  are 
parallel  to  the  x-axis. 

The  description  in  the  primed  coordinate  system  is  completed 
oy  a  consideration  of  the  thermal  stress  parameters.  In  the  fiber 
coordinate  system  the  thermal  stresses  are  given  by 


y3aT 


y^aT 


(5.12) 


wnere  aT  is  the  temperature  increment,  and  the  stress- temperature 
coupling  coefficients  are  given  in  terms  of  the  linear  coefficients 
of  thermal  expansion  by 


Yls  (CU  +  C12,ax  +  C13az 
t3=  2C13ax  +  C33az  • 


The  nonvunishing  components  of  tnis  second  rank  tensor  in  the 
primed  coordinate  system  are 


Y11  =  yl 

s  2  2 

y22  =  Ylcos  9  +  Y3sin  9 
,  .  2  .  2 
y33  =  'rlsin  9  +  ^3C0S  0 

y^3  =  Y32  —  ( Y3  -  Y|)  sinecose 


(5.14) 
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In  summary,  the  stresses  are  given  in  terms  of  the  strains  in 
tne  rotated,  or  laminate  coordinate  system  (LCS),  by 


all  =  4l£ll  +  C12e22  +  C13e33  +  2C14e23  '  y11aT 

°22  =  C12ell  +  C22e22  +  C23e33  +  2C24£23  '  y22aT 

°33  =  C13ell  +  C23e22  +  C33e33  +  2C34£23  '  y33aT 

°23  =  C14ell  +  C24e22  +  C34e33  +  2C44£23  '  Y23*T 

a13  =  2C56£13  +  2C56£12 
a12  =  2C56£13  +  2C66£12  * 


(5.15) 


wnere,  for  convenience,  the  primes  have  been  omitted.  If  we 
consider  the  stresses  as  a  six  component  column  matrix,  then  the 
above  relations  can  be  written  in  matrix  form  as 


o  =  Ce  -  y  aT  , 

and  tne  inverse  relations  giving 
stresses  are 

e  =  C  +  C  AT  , 

or  equivalently  as 

e  =  Sa  +  a  aT  , 

i  .e. , 

£11  =  Sllall  +  S12°22  +  613a33 
e22  =  SI2°11  +  S22a22  +  S23a33 
e33  =  S13all  +  S23a22  +  S33a33 
e23  =  S14a 11  +  S24°22  +  S34°33 
£13  =  2S55°13  +  2S56a12 
£12  =  2S56a13  +  2S66a12 


(5.16) 

the  strains  in  terms  of  the 


(5.18) 

2S14a23  +  allAT 
2^2 ^ +  ^22^^" 

2^34a23  +  a33A"*"  (5.19) 

2S44a23  +  “23aT 
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where  the  inverse  S  of  the  matrix  C  is  easily  obtained  using  the 
computer  (note  again  that  the  matrix  S  and  the  matrix  C  are  not 
symmetric  wnen  the  strains  are  defined  by  Eq.  (5.3)). 


5.2  THERMAL  EXPANSION  COEFFICIENTS 

Here  we  shall  De  concerned  with  the  thermal  expansion  coef¬ 
ficients  of  a  laminate  containing  equal  volume  fractions  of  plies 

with  fibers  oriented  in  one  or  the  other  of  two  directions.  The 
z-axis  of  tne  laminate  coordinate  system  is  taken  as  the  bisector  of 
the  angle  between  the  fiber  directions,  as  indicated  in  Figure  5.2. 

The  elastic  constants  in  the  LCS  of  a  ply  with  fibers  in  the 
z^  direction  are  given  directly  by  the  relations  (5.11);  and  the 

elastic  constants  of  a  ply  with  fibers  along  z^  may  be  obtained 

from  the  connections  (5.11)  by  replacing  sine  by  -  sine,  i.e.. 


-(2) 

r(l) 

°14  = 

-  C14 

-(2) 

r(l) 

-24  = 

-  b24 

~(2) 
^34  s 

r(l) 
"  U34 

,(2) 

r(l) 

^56  = 

-  C56 

(2) 

(1) 

y23  = 

'  y23 

(2) 

(1) 

a23  = 

‘  a23 

all  the  other  parameters  being  equal.  If  the  laminate  thickness 
(x-direction)  is  very  much  less  than  the  lateral  dimensions  (y,z 
directions),  then  away  from  the  edges  the  shear  stresses  between 
plies  will  vanish,  and  the  appropriate  boundary  conditions  for 
determining  tne  thermal  expansion  coefficients  are 
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Fiber 

z2 


Figure  5 


.2.  Laminate  geometry. 


7/ 


(5.21) 


(5.22) 


(5.23) 


(5.24) 


and  (equal  volume  fractions), 

(1)  +  (2)  n 
a  zz  °zz  =  0 

r.^  +  -  0 


( 1 )  +  (2) 

n  yz 


(5.25) 


(5.26) 


(5.27) 


The  conditions  (5.22),  (5.23),  and  (5.24)  give  the  three 

equations. 


(s(l)  ♦  s(2,ln(1)  ♦  (s(1)  +  c(2)^(D 
S23  S23  ,ayy  (S33  S33  )azz 

+  +  ,(2).  (1)  .  ,  (1)  (2) .  T  . 
2^S34  S34  *°yz  ^a33  “  a33  )aT  =  0 

K(lj  +  <-(2).  (1)  +  ,-(l)  +  .(2),  (1) 

)  a>/w  +  (S23  +  S23  )azz 


22  22  ,uyy 


+  ,K(1)  +  <j(2K  (1)  +  ,  (1)  (2).,  „ 

^(->24  ^24  °yz  'a2Z  ~  a22  ^  =  u 

iod;  +  c(^)^  U)  +  (cd)  +  c(2)v  (i) 

S24  S24  °yy  (S34  S34  )azz 

+  ?(Al)  +  ^(2) .  (1)  +  ,  (1)  (2)  v  T  . 

+  2(S44  S44  )ayz  ♦  (a23  -  a23  )aT  =  0 


(5.28) 


(5.29) 


(5.30) 


where  use  nas  been  made  of  (5.21),  (5.25),  (5.26),  and  (5.27).  Since 

c(2)  <■  ( 1 )  (2)  (1)  . 

34  ~  "  ^34  ’  a33  "  a33  ’  etc*» 


the  first  two  relations  (5.28)  and  (5.29)  give  the  intuitively 
obvious  result 
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(1)  (1)  (2)  (2)  n 

a  =  a  s  a  =a  =u 

yy  22  yy  zz 


and  tne  third  condition  (5.30}  gives 

aA  1)  (1)  +  o  (1)aT  n 
4S44  ayz  2a23  aT  =  0  * 


The  z-component  of  the  thermal  strain  is  just 
(1)  _  yAl)  (1)  +  (1)  aT 


ZZ  zz 


'34  yz 


‘33 


whicn  witn  tne  use  of  (5.32)  gives  the  desired  result 

c(l> 


zz  '34  (1)  +  (1) 

-r—t-  n _  ▼  r*  .  _ 


az  =  TT  =  “  7TT)  a23 
544 


‘33 


Similarly, 


?,(1)  (1)  +  (1)  aT 

eyy  =  2S24  °yz  °22  aT 


which  witn  tne  use  of  (5.32)  gives 


-d) 

>24 


(1)  +  (1)  , 


TIT  a23  a22 

*44 


and  finally,  one  obtains. 


ax  = 


,(1) 

*14  (1)  +  (1) 

'  TTU  a23  ali 

544 


(5.31) 

(5.32) 

(5.33) 

(5.34) 

(5.35) 

(5.36) 

(5.37) 


Rememoer  tnat  in  the  above  expressions  the  S . .  and  a.,  are 

■  si  *  J 

elements  of  the  tensors  in  the  laminate  coordinate  system. 

The  elements  S.  .  may  be  obtained  in  two  ways:  by  computing 

J 

the  elements  of  tne  matrix  C  and  taking  the  inverse,  or  Dy  direct 
computation  from  the  elements  in  the  fiber  coordinate  system  using 
tne  transformation  Eqs.  (5.11) 
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For  example. 


3  3 

-  S^sin  ecose  +  S-^sinecos  e 

3  3 

+  Seisin  ecose  -  sinecos  e) 

3  3 

+  2S4^(sin  ecose  -  sinecos  e)  , 


(5,38) 


wnere 


11  = 

1/Ex  -  S22 

12  = 

-  v  /£  = 

xy  x 

S21 

13  = 

-  \z/Ez  ■ 

S31 

33  = 

1/Ez 

'44  = 

1^(4C44)  = 

1/(4 

66  = 

1/(4C66)  . 

1/(4 

(5.39) 


are  tne  nonvanishing  components  in  the  FCS  in  terms  of  the  Young's 
moduli  E  E  and  the  Poisson’s  ratios  v  ,  v  (contrac- 

A  4  Aj  A4 

tion,  force).  The  quantities  G^,  Gy  are  the  standard  axial  and 
transverse  shear  moduli. 


5.3  Elastic  parameters 

Tne  elastic  (Young's)  modulus  in  the  z-direction  for  the 
laminate  is  obtained  by  applying  the  conditions 


(1) 

(2) 

zz 

=  ezz 

—  c 

(1) 

(2) 

yy 

eyy 

(i) 

-  e(2) 

yz 

yz 

(l) 

(2) 

=  0 

'xx 

'  axx 

(1) 

yy 

+  (2) 

+  a 

yy 

=  0 

(5.40) 

(5.41) 

(5.42) 

(5.43) 

(5.44) 
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=  0 


(5.45) 


where  (5.44)  ana  (5.45)  imply  equal  volume  fractions  of  the  two 
types  of  ply. 

The  condition  (5.41)  gives 


c(l)  U)  +  AD  (1)  +  oc(l)  (1) 

S22  °yy  S23  °zz  2S24  °23 

o(2)(2)  +  _(2)  (2)  +  ?.(2)  (2) 

=  b22  ayy  S23  Jzz  2S24  a23  ’ 


and  (5.42)  gives 


(5.46) 


c(l)  (1)  +  AD  (1)  +  .,(1)  (1) 

S24  ayy  S3*  °2z  2S44  a23 

s(2)  (2)  +  -(2)  (2)  +  p-(2)  (2) 

'  S24  °yy  S34  °zz  2S44  a23  ’ 


(5.47) 


wnere  use  has  been  made  of  (5.43).  With  the  use  of  (5.44)  and 
(o.45)  the  above  relations  become 


(  (i)  +  A2)  >(1)  +  ADJD  AZ)  (2) 

S22  S22  yy  S23  azz  '  S23  °zz 

+  2(S<JJ  +  S24))a23)  =  0  ’ 

(AD  +  s(2)la(1)  +  <;(i)  (1)  c(2)  (2) 

(S24  S24  ayy  S34  azz  *  S34  azz 

+0/c(l)+r(2),(l)  ~ 

2(S44  S44  ,a23  =  0  » 


(5.48) 


(5.49) 


and  (5.40)  gives  the  two  relations 

(1)  (1)  +  (1)  (1)  +  _(1)  (1) 


'23  yy 


9  .  U 

3J  zz 


34  u23 


A  2)  (1)  +  -(2)  (2)  ?.(2)  (1) 

^23  ayy  S33  azz  “  2S34  a23  =  £o 


(5.50) 


(5.51) 


The  relations  (5.48),  (5.49),  (5.50),  and  (5.51)  constitute  a  set  of 

four  equations  in  the  four  variables  o^,  o^,  a^,  and  oi1^. 

yy  zz  zz  23 
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Equation  (5.48)  simplifies  to 


>A  1)  (1)  +  AD,  (1)  (2h  n 

i%22  °yy  S23  Uzz  '  °zz  >  *  0  > 


(5.52) 


and  (5.49)  reduces  to 

eU),  (1)  +  (2h  +  4c(l)  (1) 

S34  0  zz  0 zz  4S44  a23 


(5.53) 


The  relation  (5.51)  is  also 

c(l)  (1)  +  c(l)  (2)+  ~(1)  (i) 
"  S23  °yy  S33  °zz  2S34  °23 


(5.54) 


which  when  subtracted  from  (5.50)  gives 

,c(l)  (1)  +  c(l)/  (1)  (2),  n 

2S23  °yy  S33  {azz  "  °zz  }  =  0  ’ 


(5.55) 


and  a  consideration  of  5.52)  and  (5.55)  gives  the  results, 

(1)  (2)  n 

o  =0  =  U 

yy  yy 

a(1)  -  a(2) 

zz  “  zz  * 


(5.56) 


wnicn,  it  must  be  adnitted,  seem  very  appropriate. 


The  relation  (5.53)  then  becomes 

,c(l)  (1)  +  .AD  (1)  n 

‘iS34  °zz  4S44  °23  *  0 


(5.57) 


and  (5.54)  is  now 


c(l)  U)  +  9AD  (1) 

b33  °zz  2S34  a23 


(5.58) 


The  aDove  two  relations  yield  the  desired  result, 

(1)  <-(1)2 

°zz  c(l)  b34 
— 3  s33  -  11 

0  ^44 


(5.59) 
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(5.60) 


,(1)2 

1  .(1)  ^34 

r-  =  s33  -  rnr 

z  ^44 


for  the  reciprocal  of  Young's  modulus  Ez. 

Tne  elastic  modulus  in  the  y-direction  follows  from  the 
conditions 


.U)  _ 
"yy 

(2) 

eyy  ~  eo 

(5.61) 

(l) 

ezz  = 

(2) 

ezz 

(5.62) 

(1) 

^yz 

(2) 

"yz 

(5.63) 

(i) 

axx 

(2)  _ 
axx  =  0 

(5.64) 

and  for  equal  volume  fractions. 


(1)  +  (2) 
r  t  a 


zz 


zz 


(1)  .  (2)  n 
j  +  a  =  0 

yz  yz 


(5.65) 

(5.66) 


Ooing  tne  same  algebra  as  above  for  E^  gives  the  desired  relation 


o(l)  2 

1  c(l)  *24 

T-  -hi  - 

y  ^>44 


(5.67) 


for  the  reciprocal  of  the  Young  s  modulus. 

in  a  similar  fashion,  one  obtains  for  the  elastic  modulus  in 
the  direction  perpendicul ar  to  the  laminate, 


1 

F~  = 
x 


o(l)  2 
*14 


(5.68) 
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The  corresponding  Poisson's  ratios  are  easily  obtained  from  the 
aoove  results.  For  the  case  of  an  applied  stress  in  the  z  direction. 


exx 

S13azz  +  2S14°23 

ezz  = 

S33azz  +  2S34a23 

which  wi th 

the  use  of  (5.57)  gives 

S13S44  -  S14S34 

vxz 

S33S44  '  S34 

Likewise, 

£yy  = 

S23°zz  +  2S24°23  » 

and 

S23S44  ~  S24S34 

V  = 

2 

C  ^  qt- 

‘  o4 

In  order  to  obtain  an  expression 

modulus,  we 

apply  the  conditions: 

(1) 

U) 

£yz 

eyz  “  eo 

(1) 

(2) 

ezz 

■  ezz 

(1) 

(2) 

£yy 

£yy 

(l) 

0 

(2) 

=  a  =  U 

XX 

XX 

(1) 

o 

*  .  a 

ZZ 

ZZ 

U> 

,  (2)  n 
■  c  =  U 

yy 

°yy 

(5.69) 


(5.70) 


(5.71) 


(5.72) 

in-plane  shear 


(5.73) 

(5.74) 

(5.75) 

(5.76) 

(5.77) 
(5.8/) 


The  condition  (5.74)  yields  witn  the  use  of  (5.76),  (5.77),  and 

(3.78), 
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(5.79) 


c  (1)  +  c  (1)  +  cU)/  (1)  +  (2h 

S23°yy  +  S33azz  S34  (cyz  ayz  } 


=  0 


and  (5.76)  gives 


,  Ui  +  c  (1)  +  cd)/  (1)  +  (2)v 

S22ayy  S23azz  S24  {ayz  ayz  ] 


=  0 


(5.80) 


The  superscript  is  omitted  in  the  above  expression  if  the  quantity 
has  tne  same  numerical  value  for  both  types  of  ply  (e  and  -  e  fiber 
orientation) . 

From  (5.73)  we  obtain  two  relations 


,U)  (1)  +  <-(1)  (1)  +  (1) 

S24  ayy  S34  °zz  2S44<’yz 

c(2)  (2)  .  AZ)  (2)  .  (2) 

S24  °yy  S34  °zz  2S44°yz 


=  e 


=  e 


(5.81) 

(5.82) 


Subtracting  ( 5.82)  from  (5.81)  leads  to  the  conclusion  that 

(1)  (2) 

a  =  a  , 

yz  yz 

whereas  adding  (5.82)  and  (5.81)  gives  the  useful  relation, 

Al)  (1)  +  ,.(1)  (1)  +  (1) 

S24  ayy  S34  azz  ^S44ayz  =  eo 


(5.83) 


(5.84) 


wnicn  togetner  with  (5.79)  and  (5.80)  form  a  set  of  tnree  simultane¬ 


ous  linear  equations  in  the  variables 
(1) 


(1) 

Jyy  * 


.(1), 

zz 


and 


ayz  *  These  may  be  solved  numerically  to  obtain  the  in-plane 
shear  modulus 


(1) 

r(effJ  ayz 

°44  =  7T~ 

o 


(5.85) 


5.4  NUMERICAL  RESULTS 

Calculated  results  for  seyeral  parameters  as  a  function  of  the 
lay-up  angle  ©  (Figure  5.2)  are  presented  for  a  laminate  whose  plies 
consist  of  i'J  pen.cr.1  oy  volume  graphite  fioers  in  a  2024  aluminum 
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matrix.  The  input  parameters  for  the  fibers  and  matrix,  as  given  in 

TaDle  5.1  below,  are  the  same  as  those  used  by  Hashin  and 
(2) 

Humphreys  in  a  recent  report. 

TABLE  5.1 

INPUT  ELASTIC  CONSTANTS  FOR  FIBER  AND  MATRIX 


Elastic 

Constant 

T-50 

Graphite  Fiber 

2024-T4  A1 

Matrix 

Ez(GPa) 

388.2 

72.4 

Ex(GPa) 

7.6 

72.4 

14.9 

27.2 

C66(GPa) 

2.6 

27.2 

v  ** 

xz 

0.41 

0.33 

v„  , 

0.45 

0.33 

xy  i 

<*_(  r1) 

-  0.68  x  10'6 

22.5  x  10~6 

°x(  0  J 

9.74  x  IQ"6 

22.5  x  10~6 

**Tne  subscripts  are  (contraction,  force). 

The  PRUFC  code^  was  used  to  calculate  the  elastic 
properties  of  a  ply  in  the  fiber  coordinate  system,  with  the  results 
as  given  in  Table  5.2  (the  PRUFC  code  was  subsequently  revised  as 
discussed  in  Section  6). 

The  linear  tnermal  expansion  coefficients  for  the  laminate,  as 
calculated  with  Eqs.  (5.34),  (5.36),  and  (5.37)  are  plotted  as  a 
function  of  ply  angle  in  Figure  5.3.  It  will  be  noted  that  tne 
axial  tnermal  expansion  has  a  minimum  at  about  20  deg  (40  deg  total 
angle  between  fiber  directions).  It  is  also  comforting  to  note  that 
the  calculated  results  exhibit  the  required  symmetry  about  45 
degrees. 
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The  elastic  moduli  along  the  three  laminate  axes  are  plotted 
in  Figure  5.4,  and  two  Poisson's  ratios  are  given  in  Figure  5.5  as  a 
function  of  lay-up  angle.  The  in-plane  shear  modulus  for  the 
laminate  is  shown  in  Figure  5.6. 


LINEAR  COEFFICIENT  OF  THERMAL  EXPANSION  (10" 


Figure  5. 3. Calculated  linear  thermal  expansion  coefficients  for  a 
cross-ply  metal-matrix  laminate  as  a  function  of  lay-uD 
angle. 
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6.  TRANSVERSE  ELASTIC  PROPERTIES  OF  A  UNIDIRECTIONALLY 
REINFORCED  COMPOSITE  AND  THE  PRUFC  CODE 


r  2'| 

In  a  recent  report,  J  Hashin  and  Humphreys  present 

numerical  results  for  a  typical  graphite/aluminum  composite  as 

calculated  using  a  finite-element  code  (ANSYSJ  to  determine  tne 

stress-strain  state  in  the  unit  cell  of  a  model  in  which  the  fibers 

are  arranged  in  a  perfect  hexagonal  array.  The  advantage  of  their 

method  is  tnat  it  is  expected  to  give  a  very  accurate  representation 

of  tne  behavior  of  a  composite  in  terms  of  the  properties  of  the 

fiber  and  matrix  components.  The  disadvantage  is,  as  stated  in  the 

report,  "The  ANSYS  computer  code,  ....  is  very  expensive  to  run  for 

sucn  cases."  Consequently,  it  is  of  interest  to  compare  these 

f31 

finite  element  results  with  those  obtained  with  PRUFCL  and  the 
S-CUBED  plastic-flow  code,^  which  are  inexpensive  to  run  but  use 
the  concentric-cylinder  approximation. 

He  present  first  a  comparison  of  calculated  elastic  properties 
of  a  graphite/aluminum  composite  as  obtained  from  the  properties  of 
the  constituents.  In  the  example  given  in  Ref.  2,  the  properties  of 
the  fiber  and  aluminum  matrix  are  as  given  in  Table  6.1. 

A  comparison  of  tne  composite  parameters  as  calculated  with 
PRUFC  and  those  calculated  by  Hashin  and  Humphreys  for  a  30  percent 
by  volume  of  fioers  is  given  in  Table  6.2  below.  Also  given  are 
values  obtained  from  the  "Composite  Cylinder  Assemblage"  (CCA)  model 
of  a  composite  (Hashin,  Ref.  8). 

Tne  only  significant  differences  are  in  the  transverse  proper¬ 
ties  £y.  Gyp  and  v-p  (only  two  of  which  are  independent,  be¬ 
cause  of  tne  transverse  isotropy).  Not  all  the  transverse  pro¬ 
perties  differ,  however.  If  one  considers  the  transverse  bulk 
modulus  K  wnose  reciprocal  is  given  by 
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TABLE  6.1 


CONSTITUENT  PROPERTIES  FOR  A  GRAPHITE/ALUI1INUM  COMPOSITE 
Elastic  T-50  Graphite  2Q24-T4  A1 


Constant 

Fiber 

£A(GPa) 

388.2 

ET(GPa) 

7.6 

GAT(GPa) 

14.9 

Gyy(GPa) 

2.6 

\J  * 

AT 

0.41 

VTT  , 

"A(*C  1 

0.45 

-0.68 

^(V1) 

9.74 

A  =  Axial  (Longitudinal) 
T  =  Transverse 


Matrix 

72.4 

72.4 

27.2 

27.2 


0.33 

0.33 

x  10"6 

22.5  x 

10'6 

LO 

1 

o 

X 

22.5  x 

io~6 

*Tnis  Poisson's  ratio  is  tne  transverse  contraction  fur  an  axial 
stress. 


TABLE  6.2 

CALCULATED  PROPERTIES  FOR  A  GRAPHITE/ALUMINUM  COMPOSITE 
(30  Percent  by  Volume  of  Fibers) 


Elastic 

Constant 

PRUFC 

Finite-Element 

Method 

CCA 

EA(GPa) 

167.18 

167.14* 

-- 

Er(GPa) 

33.06 

42.26 

41.78 

Jyy(GPa) 

10.86 

15.13 

14.99 

22.82 

23.20 

22.87 

VTT 

VAT  -1 
aA(  C  ) 

aT(°C_1) 


0.522 
0.338 
6.29  x 


10 


25.65  x  10 


-6 

-6 


0.39o 
0.340 
6.36  x 


10 


-6 


25.69  x  10' 


0.394 
0.338 
6.36  x  10 
25.65  x  10 


-6 

-6 
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the  values  as  obtained  by  the  PRUFC  code  and  the  finite  element 
method  are,  respectively, 

k  =  38.19  GPa  (PRUFC) 
k  =  38.73  GPa  (F.E. )  , 

a  good  agreement.  The  difficulty  is  that  the  concentric-cylinder 
approximation  to  a  composite  does  not  provide  a  unique  definition  of 
tne  transverse  shear  modulus. 

In  the  PRUFC  code  as  documented  in  Ref.  3,  the  transverse 
shear  modulus  is  obtained  from  a  consideration  of  the  deformation  of 
concentric  cylinders  subjected  to  a  stress  field  at  its  outer 
surface  the  same  as  that  acting  on  a  fictitious  cylindrical  surface 
in  a  nomogeneous  material  subjected  to  tne  stress  field 


a 


XX 


ayy  =  °o  _  cons^* 


e 


Z 


=  0 


(see  Appendix  B  of  Reference  9).  For  the  fictitious,  homogeneous 
cylinder,  of  course,  the  resulting  strain  will  be  uniform  throughout 
its  volume  and  is  related  to  the  stress  field  by 

°J  =  2iTTexx  •  (6'21 

involving  only  tne  transverse  shear  modulus,  Gyy.  For  the  concen¬ 
tric  cylinder,  however,  for  which  an  exact  solution  of  the  equations 
of  elasticity  is  available,^  the  effective  strain  6x/x  (see 
Figure  6.1)  at  the  outer  surface  will,  in  general,  no  longer  be 
constant  DUt  will  vary  witn  tne  angle  e.  In  fact,  for  some  com¬ 
binations  of  properties  and  fiber  volume  fractions,  the  effective 
strain  6 x/x  may  go  from  positive  to  negative  as  e  ->  90  degrees. 
Thus,  if  the  relation  (6.2)  above  is  used  to  obtain  the  shear 
modulus  from  tne  concentric-cyl inder  model,  the  value  of  e  used 
must  be  some  sort  of  average.  The  same  problem  of  nonuniqueness 
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r  q  i 

occurs  with  Hasnin's  CCA  model,  J  the  value  for  the  transverse 
snear  modulus  quoted  in  Table  6.2  being  the  greater  of  calculated 
upper  and  lower  bounds  (see  page  295  of  Reference  8). 

In  the  original  version  of  PRUFC,  a  volume  weighted  average 
was  used  for  exx  in  Eq.  (6.2)  to  obtain  the  transverse  shear 
modulus,  i.e.. 


(6.3) 


ax  being  tne  displacement  at  the  outer  surface  r  =  r2.  This 
procedure  results  in  the  values  for  the  transverse  properties  given 
in  Taole  6.2.  Since  the  transverse  shear  modulus  obtained  in  this 
manner  is  significantly  below  that  as  calculated  by  the  finite- 
element  metnod,  which  is  presumably  the  better  value,  several  other 
averaging  alternatives  have  been  examined.  Results  are  given  in 
Figure  6.2,  where  the  shear  modulus  is  plotted  versus  fiber  volume 
fraction. 

The  use  of 


(6.4) 


results  in  a  decrease  in  tne  calculated  shear  modulus  relative  to 
the  volume  weighted  average  (6.3),  and  a  greater  discrepancy  with 
riashin  and  Humphrey's  result  for  the  30  percent  composite.  A  closer 
agreement  with  the  finite  element  result  is  obtained  if  one  uses  the 
linear  average 
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i  i  i  i  i  I  i  i  r 


0  0.5  1.0 


FIBER  VOLUME  FRACTION 

Figure  6.2.  PRUFC  results  for  the  transverse  shear  modulus  of  a 

graphite/aluminum  composite  as  a  function  of  fiber  volume 
fraction  for  various  possible  averaging  procedures  using 
stress  boundary  conditions.  Also  showp  is  the  finite- 
element  result  of  Hashin  and  Humphreys^]  for  30  percent 
by  volume  of  fibers. 


but  the  PRUFC  calculation  is  still  low.  Finally,  the  use  of 


-(H 


(6.6) 


in  Eq. (6.2)  results  in  a  PRUFC  shear  modulus  greater  than  that 
obtained  with  the  finite  element  method.  The  use  of  Ux/x)  at  45° 
suggests  itself  because  tnis  is  also  the  value  of 


(?-£)  • 


wnicn  turns  out  to  be  independent  of  e.  Calculated  results  for 
v  and  £  using  the  various  values  of  l  are  given  in 

Ajr  a  XX 

Figures  6.3  and  6.4. 

An  alternative  approach  is  to  specify  the  displacement  on  the 
outer  surface  of  the  concentric  cylinder  model  and  to  solve  the 
equations  of  elasticity  for  the  resulting  stress  field.  As  outlined 
in  Ref.  9,  the  radial  and  hoop  displacements  are  of  the  form 


U^a,(r,e)  =  u^(r)  cos2e 
V^a)(r,e)  =  v^ai(r)  sin2e  , 


(6.7) 


(6.8) 


and  the  condition  now  that  the  displacements  on  the  outer  surface. 


6X  ay 

~  =  "  y  =  e0 


constant 


(6.9) 


be  the  same  as  tnat  on  a  cylindrical  surface  within  a  homogeneous 

material  subjected  to  the  uniform  strain  field  e  =  -e  = 

xx  yy 

Eq,  is  satisfied  by  requiring  that 


v(2)(r2) 


u(2)(r2)  . 


(6.10) 


The  stress  boundary  conditions  (B.9)  and  (B.10)  of  Ref.  9  are 
replaced  by 
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I  I  I  I  I  I 


•  Finite  element  result 

©  Sx  s  (“X/XW  Eq'(6-4) 

(T)  Volume-weighted  average,  Eq.  (6.3) 
(T)  Linear  Average,  Eq.  (6.5) 

0  =xx  =  <5x/x>9=45°’  Eq'  (6>6) 


FIBER  VOLUME  FRACTION 

Figure  6.3.  PRUFC  resuits  for  the  transverse  Poisson's  ratio  of  a  graphite/ 
aluminum  composite  as  a  function  of  fiber  volume  fraction  for 
various  possible  averaging  procedures  using  stress  boundary 
conditions.  ALso  shown  is  the  finite-element  result  of  Hashin 
and  HumphreysL?]  for  30  percent  by  volume  of  fibers. 


US 


(GPa) 


FIBER  VOLUME  FRACTION 

Figure  6.4.  PRUFC  results  for  the  transverse  elastic  modulus  of  a  graphite/ 
aluminum  composite  as  a  function  of  fiber  volume  fraction  for 
various  possible  averaging  procedures  using  stress  boundary 
conditions.  ALso  shown  is  the  finite-element  result  of  Hashin 
and  Humphreysl^J  for  30  percent  by  volume  of  fibers. 


1 


(6.11) 


J2),  r  v 

u  (r2'  =  eo  * 
v  ( r2  ^  =  ~  £o  * 


(6.12) 


the  other  conditions,  (6.5)  -  (B.8)  remaining  the  same.  Equations 
(6.11)  and  (6.12)  follow  from  the  relations  for  the  cartesian 
components  of  displacement  in  terms  of  tne  radial  and  hoop 
components. 


<sx  =  u(r)  cos2e  cose  -  v(r)  sin2e  sine  (6.13) 

ay  =  u( r)  cos2e  sine  +  v(r)  sin2e  cose  .  (6.14) 


As  before,  the  boundary  conditions  lead  to  a  set  of  six  linear 
equations  for  the  coefficients  of  the  elastic  solutions,  Eqs.(B.27) 
of  Ref.  9,  the  bottom  two  rows  of  the  matrix  K  being  replaced  by 


K51  = 

1. 

II 

CM 

LO 

1,  K 

53 

=  1, 

K54  = 

1, 

K55  = 

0,  K 

56 

=  o  , 

K61  = 

-1. 

K62 

1 

*  '  1 

(3 

n 

♦n(2>) 

m —  • 

II 

VO 

1 

'  2 

(1  - 

U). 
n  ) , 

K 

64  =  1  » 

K65  * 

o. 

K6o  = 

0  , 

where  the  applied  strain  and  outer  radius  have  been  normalized  to 
e  *  1,  =  1.  The  coefficient  Kg3  as  given  in  Ref.  9  is  in 

error.  It  should  be 


2K1 
n  )-=- 
rl 


From  £q.  (6.13)  of  Ref.  9,  the  radial  stress  at  the  surface  of  the 
outer  snel 1  i s  gi ven  by 
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(1,9) 


and  from  (B.15) 


=  Cc[^(A2  ♦  3B2  -  C2  -  3D2)  -  cj2}J  cos  2e 

5  cos2e 
rro 

the  shear  stress  is 


=  C 


(2) 

66 


A2  '  I 


3  (3  +  n(2)) 


n 

T Tf 


*  j  (1  -  nl2,)C2  -  3D2  -  1 


b2 

sin2e 


=  sin2e 
reo 


(6.16) 


(6.17) 


(Actually,  the  aoove  quantities  are  stresses  per  unit  of  applied 
strain,  a/z^).  As  in  Eq.  (6.2)  above,  the  transverse  shear 
modulus  will  be  oDtained  from-  the  relation 


"xx  *  ^TT‘0 


(6.18) 


wnere  ?xx  is  some  sort  of  average  at  the  outer  surface.  In 
arriving  at  an  average,  we  consider  tne  x-component  of  the  surface 
traction  rather  tnan  °xx(r,9)  itself.  In  general,  the  force  per 
unit  area  on  a  surface  with  normal  unit  vector  n  is  given  by 


t 


and  for  tne  cylindrical 


surface  n  =  r  so  that 


A  A 

ra  +  Qa 
rr  re 


and  tne  x-component  is 

fx  =  arrcose  -  are  sine 


(6.19) 


(6.20) 


(6.21) 


(Note  tnat  a  itself  also  depends  upon  a  ).  The  stress  per 

A  A  OO 

projected  area  perpendicular  to  tne  x-axis  is 
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(6.22) 


0  =  i  a  -  i  a  (6.24) 

xx  3  rro  3  reo  v  ' 

where  a  „  and  <j  „  are  defined  by  (6.16)  and  (6.17). 
rro  reo 

A  comparison  of  the  calculated  transverse  shear  modulus  using 
the  linear-average,  displacement  boundary  conditions  with  that  from 
the  linear-average,  stress  boundary  conditions  is  given  in  Figure 
6.5.  Also  plotted  is  the  arithmetic  mean  of  the  two  shear  moduli, 
wnicn  altnougn  still  hign  witn  respect  to  the  finite  element  result 
appears  to  be  the  most  plausible  way  to  obtain  the  transverse 
properties  using  tne  concentric-cylinder  approximation,  and  the 
PRUFC  code  has  been  revised  accordingly.  A  comparison  of  results 
for  the  30  percent  composite  using  the  revised  PRUFC  code  witn  those 
obtained  by  the  finite-element  method  is  given  in  Table  6.3. 

F21 

Hashin  and  Humphreys  also  present  calculated  results  for 
a  thermal  process  in  which  the  yield  strength  of  the  aluminum  matrix 
material  in  a  30  percent  composite  is  exceeded.  The  yield-strength 
model  they  used  is  given  as  Figure  6.6,  which  is  reproduced  from 
Ref.  2.  This  isotropic  hardening  model  was  used  in  the  S-CUBEO 
residual  stress  code^  along  with  the  fiber  and  matrix  properties 
as  given  in  Table  6.1,  and  two  of  the  calculations  were  repeated. 
Figure  6.7  gives  the  axial  strain  as  calculated  with  the  S-CUBED 
model^i_l  for  a  30  percent  composite,  initially  stress  free  at 
37i°C,  as  it  is  cooled  to  room  temperature  and  reheated.  Using  the 
isotropic  hardening  model  of  Figure  6.6,  tne  matrix  first  reaches 
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(GPa) 


Fiaure  6 . 5 . Compari son  of  the  transverse  shear  modulus  for  a  grapnite/ 
aluminum  composite  as  calculated  with  disolacement  boundary 
conditions  and  stress  boundary  conditions.  Also  shown  is  the 
finite  element  result  of  Hashin  and  Humphreys. T2] 


the  yield  surface  at  ”142°C.  Upon  reheating  from  room  temperature, 
the  matrix  immediately  falls  off  the  yield  surface  and  does  not 

TABLE  6.3 

CALCULATED  PROPERTIES  USING  THE  REVISED  PRUFC*  CODE 
FOR  A  GRAPHITE/ALUMINUM  COMPOSITE 
(30  percent  by  volume  of  fibers). 


Elastic 

Constant 

PRUFC* 

Finite-Element 

Method 

EA(GPa) 

167.18 

167.14 

ET(GPaj 

45.22 

42.26 

G-j--j-(GPa) 

16.80 

15.13 

GA-j.(GPa) 

22.82 

23.20 

VTT 

0.346 

0.396 

VAT  i 

aA(  °C  ) 

vr‘i 

0.338 

6.29  x 

25.65  x 

io-6 

10-6 

0.340 

6.36  x  10 

25.69  x  10 

*Revised  as  per  above. 

reacn  it  ayain  until  about  350°C.  The  initial  segment  from  371°C  to 

142°C  is  shown  as  dashed  because  the  present  version  of  the  S-CUBED 

code  does  not  produce  output  for  this  interval.  The  results  for  the 

T2  1 

same  process  obtained  by  Hashin  and  Humphrey sL  using  the  finite 
element  method  and  ANSYS  code  are  given  in  Figure  6.8.  The  agree- 
menc  between  the  two  types  of  calculation  appears  to  be  quite 
satisfactory. 

Calculated  results  for  the  transverse  strain  for  the  same 
tnennal  cycle  are  given  in  Figures  6.9  and  6.10.  Again  the 
agreement  is  quite  satisfactory .  We  conclude  tnat,  when  their  use 
is  appropriate,  the  revised  S-CUBED  PRUFC  code  and  concentric-cylin¬ 
der  residual -stress  model ^  will  give  results  for  metal -matrix 
composites  which  are  as  satisfactory  as  those  obtained  with  the  more 
complex  finite  element  metnod. 
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no 


Figure  6.9.  Transverse  strain  as  calculated  with  the  S-CUBED  residual  st 
for  a  metal -matrix  composite  cooled  from  an  initial  stress-f 
at  371°C.  Results  for  the  same  process  using  a  finite  eleme 
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7.  THERMAL  DEFORMATIONS  AND  STRESSES 
IN  AN  ELASTIC  PLATE  UNDER  AN 
ARBITRARY  TEMPERATURE  DISTRIBUTION 


Tne  material  in  this  section  has  previously  been  submitted  as 
a  Progress  Report  to  the  Defense  Supply  Service  (Ref.  10).  It  is 
reproduced  here  since  it  might  prove  useful  in  computing  the  per¬ 
manent  bending  deformations  of  a  metal -matrix  composite  under  X-ray 
energy  deposition. 

Tne  method  of  separation  of  variaoles  is  used  to  obtain  the 
elastic  displacements  in  a  two-dimensional  plate  subject  to  an 
arbitrary  transverse  temperature  distribution.  The  geometry  is 
illustrated  below:  x 


Tne  plate,  of  widtn  2xQ  and  length  2 zQ,  is  assumed  to  be  aniso- 
tropic-transversely  isotropic  with  the  xy  plane  being  the  plane  of 
isotropy.  The  stresses  are  then  given  by 


3U.  3Uo  3U- 

°xx  =  C11  ~  +  C12  ay"  +  C13  a T 


a  =  C 

yy 

a_  =  C 


-  yJU) 


3U. 

- *.  + 

ax 

au , 

au. 

12 

cn 

_ i  + 

ay 

C13 

o 

az 

Yj^Ttx) 

3U 

a  y 

3U0 

13 

_ L  + 

ax 

C13 

C33 

•i 

3 Z 

Y3T(X) 

(7.1) 


°zx  =  C44 


\  3Z  3 X  / 
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wiiere  u^,  u2,  u3  are  the  displacements  in  the  x,  y,  and  z 

directions,  T ( x )  is  the  applied  temperature  distribution,  and 
y3  are  the  stress-temperature  coupling  coefficients,  which  are 
given  in  terms  of  the  linear  expansion  coefficients  a^,  S3  by 

Y1  =  (C11  +  C12}  31  +  C1363 

(7.2) 

y3  =  2C13el  +  C33S3 

With  appropriate  cnanyes,  as  indicated  below,  the  method  pre¬ 
sented  here  could  be  adapted  to  an  isotropic  plate. 

In  the  case  of  plane  strain,  =  0,  tne  two-dimensional 
problem  of  interest  here  is 


3Un 

3U-, 

axx 

=  Cu 

_ L  + 

3X 

C13 

J 

3  Z 

■  rL  T(x) 

3U. 

3U0 

azz 

=  C13 

_ _  + 

3X 

C33 

j 

Jz~  ‘ 

•  y3  T(x) 

azx 

=  C44 

(3Ui 
\  3  Z 

3X  /  * 

The  equations  for  elastic  equilibrium. 
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3  o  3a 

XX  x  XZ 

~3X~  “3 T  =  0 
3Z  3X 


become  with  the  use  of  (7.3) 


2  2 
'44  gz2  L13  L44  3X3Z  =  Y1  3x 


'44  *  <C13  *  C44>  '3HI  *  0  ' 


(7.5) 


(7.6) 


In  tne  next  two  sections,  we  outline  the  solution  of  the  set 

of  equations  (7.6)  subject  to  the  condition  that  the  shear  and 

nonnal  stresses  vanish  at  the  surface  of  tne  plate.  The  formal 

solution  itself  is  exact;  the  numerical  evaluation,  however,  is 

approximate  in  the  sense  tnat  the  Fourier  coefficients  are  obtained 
from  the  solution  of  a  system  of  coupled  equations  of  the  form 


Clm 

E 

bmn 

Aln 

Al„ 

-  fn 

E 

m 

^nm 

Clm  ■ 

> 

and 

Clm 

are 

the  desired 

n’ 

bmn’ 

^nm 

are 

constants. 

(7.7) 


erations  limit  one  to  the  use  of  a  few  Hundred  terms. 

In  Section  7.1  we  outline  the  solution  for  an  applied  tempera¬ 
ture  distribution  which  is  an  odd  function  of  x,  and  in  Section  7.2 
the  solution  appropriate  for  an  even  function.  The  displacements 
corresponding  to  an  arbitrary  temperature  function  are  obtained  from 
a  superposition  of  the  two  solutions. 
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Numerical  results  are  presented  in  Section  7.3.  In  parti¬ 
cular,  we  present  calculated  shear  strains  for  a  discontinuous 
applied  temperature  distribution  of  the  fora 

(°  .  -*o  <  *  <  °-8  Xo 

T(x)  =< 

,  0.8xo  <x<x0 

Tnis  problem  is  mathematically  equivalent  to  the  uniform  heating  of 
two  Donded  plates  which  have  the  same  elastic  constants  but 
different  coefficients  of  thermal  expansion. 


7.1  ODD  TEMPERATURE  FUNCTION 

This  part  of  the  solution  has  been  reported  in  Reference 
9,  but  in  order  that  this  report  be  self-contained,  is  repeated 

nere.  For  an  odd  temperature  distribution,  the  appropriate  combina¬ 

tion  of  functions  to  represent  the  elastic  displacements  is 

u^x.z)  =  ^2  Bn(x)  cos  kz  ♦  ^2  Dm(z)  cos  J’x 

A  “  (7.9) 

Uj(x,z)  =  An(*)  sin  kz  +  2-r  cm(z)  sin  ^x 

n  m 


where 


x  =  (n  +  j)  ,  n  =  0,  1,  2,  . .. 
4o 

j  =  (in  ♦  j)  j-  ,  m  =  0,  1,  2,  . .. 


(7.10) 


and  An(x)  and  Cm(z)  are  odd  functions  of  their  respective  vari- 
aoles,  and  Bn(x)  and  Dm(z)  are  even  functions.  The  functions 
An(x),  Bn(x),  Cm(z)  and  Dffl(z)  are  subject  to  the  following 
constraints  as  determined  from  the  boundary  conditions: 

1)  Zero  shear  stress  along  ( xQ , z) 

An‘*o>  -  *W  *  0  17-U) 
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2)  Zero  shear  stress  along  (x,zQ) 


jw  *  d;<zo>  -  -  X?  r  (-1  in /°ca;(x) 

on  •'rv 


-  kB  (x) ]  cosjxdx 
n 


3)  Zero  normal  stress  along  (xQ,z) 


(7.12) 


*  chb;uo>  ■  -  r  ?  |-1)”  J  [C13C;(Z> 
o  m  J  q 

-  jciiDm(z}]  cos  kzdz  +  ~  ^  T(xq)  (7.13) 


4)  Zero  normal  stress  along  (x,zQ) 


C,X'(zn)  -  C.JDjz 
Jj  m  o  lo  m 


°|s^3  t 

0  Jn 


T(x)  sin  jxdx  (7.14) 


The  differential  equations  for  An(x),  Bn(x),  Cmlz), 
Dm( zJ  are  obtained  from  a  substitution  of  equations  (7.9)  into 
equations  (7.6).  The  resulting  equations  for  An(x)  and  Bn ( x )  are 


C^Al'tx)  -  C33kzAn(x)  -  (Ct,  +  CAA)  kB'(x)  =  0 


'13  44'  “V 


CllBn,(x)  *  C44k  Bn(x)  +  (C13  +  C44>  kAn(x)  =  0 
and  for  Cm(z)  and  Dm(z)  are 

C33Cm ‘  ^ z >  ‘  C44j  Cm*zl  ‘  ^13  *  C44*  ^Da/zl  =  0 
C44D"(z)  -  Cl/°.lz>  *  IC13  *  C44>  KU) 


.i,  Y - 

X„  ’1  I  3X 
0  Jn 


cos  jxdx  . 


(7.15) 


(7.16) 
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The  appropriate  solutions  are 


sinhU^kx)  sinh(a2kx) 

A(j(x)  =  Aln  cosn(a^x^)  +  A2n  cosh(a2kxo) 


coshU^kx)  cosh(a2kx) 

Bn(x)  -  C0Sh(a^kx0)  +  ^2n^2  cosh(a2kxQ) 


(7.17) 


Cm(z> 


D™(2> 


sinn(bjjzj  sinh(b2jz) 

Clm  coshTo^Tz^T  +  C2m  cosh(b2jzQ 


coshfbjjz)  coshfb^z) 

ClmRl  cosh(b1jzQ)  +  C2mR2  cosiUbgj’zJ) 


i  ti  y 

~°h7  1 


3T  .  . 

—  cos  jxdx  , 


(7.18) 


where  Aln,  A2n,  Clm,  C2m  are  constant  expansion  coefficients 
to  be  determined  from  the  Doundary  conditions  (7.11  -  / .  14) .  The 
hyperbolic  cosines  in  the  denominators  were  factored  out  to  facili¬ 
tate  tne  avoidance  of  overflow  during  numerical  evaluation.  The 
quantities  a^,  a2  are  the  positive  roots  of  the  equation 


^11^44d  “  "  2C i  Ajl  -  C17)a  +  =  U 


'ir33  13  44  '  13 


33  44 


(7.19) 


_  [C44a1  -  C3^ 

1  '  un  *  "441  ai  ’  • 


and  b^,  b2  are  the  positive  roots  of 


C33C44b  *  (C11C33  "  2C13C44  '  C13)b2  +  C11C44  =  0 


(7.20) 
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* 


^33bi  ’  °44* 

«,  -  (C  t  ;  )K  •  '  -  !-2 

1  U13  U44JDi 


rtote  that  in  the  limit  of  isotropy,  equation  (7.19)  has  only  one 

2 

positive  root,  a  =  1,  and  the  solution  for  An(x)  becomes  a 
linear  combination  of  sinh(kx)  and  x  cosh(kx). 

Substitution  of  (7.17)  and  (7.18)  into  the  boundary  conditions 
(7.11  -  7.14)  yields  the  following  equations  for  the  expansion 
coefficients, 

1)  Zero  shear  stress  along  (xQ,z) 

Aln(dl'Sl)  +  A2n(a2“S2}  =  0  (7-21) 

2}  Zero  shear  stress  along  (x,zQ) 

Ciin(l  +  R^b^)  tanh  (bjJZg)  +  C2m^  +  R2b2^  tanh  ^ b2^ Zo ^ 

(7.22) 

■  -  /  (-!>"  E  k  -rkr  *  *2„  ~ir~TTr 

xo  n  [  j  +  a J YT  ^  r+apr 

3)  Zero  normal  stress  along  (xQ,z) 


Aln(C13  +  Cllslal)  tanh  {alkxo) 

+  A2n  ^C13  +  CllS2a2^  tanb  ^a2kxo^ 


Y-  (-l)n  Yj  (-1)m  J*  C 


ci3brcnRi 


x  +bjj' 


(7.23) 


C13b2“CllR2 


r  ♦  b|j‘ 


In  arriving  at  (7.23)  from  (7.13),  use  has  been  made  of  the  fact 
that  T(x)  is  an  odd  function  and  that 


x 

T(xq)  k  J]  (-l)m  j  T(x)  sin  jxdx  . 
o  m  Jn 
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4)  Zero  normal  stress  along  (x,zQ) 

Clm  (C33°1  '  C13R1)  +  C2m(C33b2  "  C13R2} 


(7.24) 


/  ^3  "  T1  j)  \  T{X)  Sin  jxdx 


prom  (7.24), 


C2n>  -  -  Q2Clm  *  Q1  j/  T(x)  s1njxdx 


(7.25) 


where 


^r- 


b  -  n  eg) 

(C33b2  "  C13R2) 


(7.26) 


n  =  C33J1  "  213R1 
2  C33b2  “  C13R2 


(7.27) 


Note  that  aoove  is  defined  differently  tnan  in  Reference  9. 
From  (7.21) 


"  Q3Aln 
dl  "  S1 


d2  -  S2 


(7.28) 


(7.29) 


Suostitution  of  (7.25)  and  (7.28)  into  (7.22)  and  (7.23)  gives  two 
sets  of  coupled  equations  for  th^  coefficients  A^n  and  C^, 


Cl,nQ5m 


j| j  T ( x )  sinjxdx  (1  +  R2D2 ^  tanb  ^b2Jzo^ 


-  ~  (-l)'n 
0 


(7.30) 

di  —  S-i  dp  -  Sp 

(-1)  k  —7 - Tj—rr  -  Qt  —7 - 7—  -  A, 


j  +atk‘ 


j  +a^k" 
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Aln^4n  =  "  z“  {_1)nQl(C13b2"CllR2 


R2^ 

m 


(-i) 


m 


^°T(x)  sinjxdx 

- 2 — 2“2 - 

k +b‘r 


^  (-l)n  (-1)m  J 


-  ^ 


m 

C13D2'C11R2 

- 2 — 2“2 — 

k*+b|r 


C13brCllRl 


'lm  * 


(7.31) 


where 

=  { 1+R1b1 )  tanh  ( b1J zq j  -  Q2(l+R2b2)  tanh  (b2jzQ)  ,  (7.32) 

Q4n  H  ^13+CllSlal)  tanh  <»ikx0>  '  V C13+CllS2a2)  tanh  ^a2kx0^ 

(7.33) 

The  coupled  sets  of  equations  (7.30)  and  (7.31)  are  of  the  form 


Clm  am  ^  bmnAln 


Aln=fn+  X  9nmClm 
m 


(7.34) 


wnicn  can  oe  solved  witn  the  following  iteration  procedure: 


AO) 

ulm 

am 

AU) 

Mln 

*  fn  * 

E 

in 

9  ci0) 

3nm  lm 

Llin 

■  dm  ' 

E 

n 

b  A(1) 
DmnAln 

AU} 

Aln 

■  fn  * 

E 

m 

9  C,(1) 
3nm  lm 
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7.2  EVEN  TEMPERATURE  FUNCTION 

For  an  even  temperature  distribution,  the  appropriate  combina¬ 
tion  of  functions  for  the  displacements  ts 

u^x.z)  =  ~y  ^  8*(x)  C0SKZ  +  ^  ^  D*(z)  sinjx 


u  j(x,z)=X  A* 


(7.36) 


A*(x)  sinkz  + 


C*(z)  cosjx  , 


where  as  before. 


k  =  (n  +  i)  ,  n  =  0,  I,  2, 
4  zo 

j  =  (m  +  i)  -J-  ,  m  =  0,  1,  2, 
o 


amd  B*n ( x )  and  C*ffl(z)  are  odd  functions  of  their  respective 
variable;,  and  A*n(x)  and  D*m(z)  are  even  functions.  The 
ooundary  conditions  impose  the  following  constraints: 


1)  Zero  shear  stress  along  (xQ,z) 


AJ'Uo'  -  k3Slxo>  -  b0  S  l-l»m  J[KU)  -  D*;U>]s1nl«dk 


(7.36) 


2)  Zero  shear  stress  along  (x,zQ) 


x 

J;;(zo)  '  °m  (zo>  3  T~  S  (-1)n  /W(xl  -  kB*U)lsinjxdx 

o  n  Jn*-  J 


3)  Zero  normal  stress  along  (xQ,z) 

cukA;<V  *cuBr(xo»  --r^i  T'V1^ 


(7.37) 


(7.38) 
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4)  Zero  normal  stress  along  (x,zQ) 


*  C13J  DS<zo>  -  3T  t3  f  Tlxl  cosjxdx  (7-J9) 

0  •'n 


The  equations  (7.6)  for  elastic  equilibrium  are  unchanged,  as 
are  tne  equations  (7.15)  for  A*n(x)  and  B*n(x).  The  equations 
for  C*ffl(z)  and  D*m(z)  become 

C33C»"lz)  -  Caa32c£(z>  *  <C,,  *  C„)  JQ^'Cz)  .  J 


'44J  m 

.2 


13  *44'  d  m 


^44^  '  ^  ~  ^iiJ  ^*(z)  "  ^n+ 


'11J  m 


13  44'  J  m 


2 

x7  Y1 


3T 


f  ^SinjXdX  * 


The  functions  A*(x)  and  B*(x)  are  now 
n  n 

cosh  (d.kx) 


cosh  (a0kx) 


u  w  jii  V  u  2  r\/\  /  VV4II  \»2IVyN/ 

A*(x)  =  A|n  cos),  (aikxQ)  +  A2n  cosh  (a2kxQ) 

sinh  (a^kx)  sinh  (a2kx) 

B*(x)  =  A|nSj_  cosh  UjkxQ)  +  A2nS2  cosh  (a2kxQ) 


and  C*(z)  and  0*(z)  become 
m  m 

sinh  (bjjz)  sinh  (b2jz) 

=  Clm  cosh  (b1JzQ)  +  C2m  cosh  (b2jzQ) 

cosh  (b.jz)  cosh  (b2jz) 

Drr/Z^  =  ClmRl  cosh  ( b^ j zQ )  +  C2mR2  cosh  (o2jzQ) 


(7.40) 


(7.41) 
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The  quantities  a^,  a 2  are  the  positive  roots  of 
CnC33a4  '  (C11C33  '  aC13C44  "  C13)  3  *  C33C44 


.  (C33ai  '  C33>  ,  j  , 

Si  -  (Cw  *  C44)a,  •  ’  -  '• 


(7.42) 


tne  same  as  for  tne  odd  temperature  function;  and  b^,  62  are  the 
positive  roots  of 

C33C44b4  "  (C11C33  "  2C13C44  "  C13,b2  +  C11C44  =  0  ’ 

but  now 


Ri  =  - 


(C33bi  "  C44) 
(C13  +  C44)bi 


,  i  =  1,  2  . 


(7.43) 


Tne  boundary  conditions  give  the  following  relations  for  the 
expansion  coefficients. 


1)  Zero  shear  stress  along  (xQ,z) 


AIn(al'V  tanh  (aikx0)  +  A2n(a2'S2)  tanh  (a2kxoJ 


(7.44) 


2  (-1)"  V  ,  ,,m(2  c.  .  1-3 

-  — -T~  i->  (‘li  J  TZZJT2 - Clm  777; 

0  m  k  +b1o  k  +b2J 


b.d-Rpl)  b2(l-R*b2) 


2)  Zero  snear  stress  along  (x,z0) 


Clm(1"Rlbl)  tanh  {bljz0*  +  C2m(1-R2b2)  tanh  (b2jzo) 


=  2  (-l)m  ^2 

3  n 


(7.45) 


.1  r  a,  (ai-S, )  a«(a^-So) 


J  +a1k 


j  +a2k‘ 
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3)  Zero  normal  stress  along  (xQz) 

Aln(C13  +  CllSlal)  +  A2n(C13  +  CllS2a2) 


n 


-  —  y  T(x  )  1 

zo  Yl  ‘  1? 


4)  Zero  normal  stress  along  (x,zQ) 
Clm{C33bl+C13Rl)  +  C2m(C33b2+C13R2) 


(7.46) 


(7.47) 


* \~i  J T(x>  cos  Jxdx  *  fo «  sinJxdx 


From  (7.46), 


AJ„  =  -  QfAJn  *  Q*T(x0) 


wi  tn 


Qf  * 


(C13*  CllVl) 

'ViM' 

zo(CU*tnS2a2l 


and  from  (7.47) , 


r* 

2m 


-  ♦  n*  i 

g3°lm  g4  J 


u 

J"  T(x)  cosjxdx 
0 


1  0 

+  ^5  7?  f 
J 


ll 

3X 


sinjxdx 


(7.48) 


(7.49) 

(7.50) 


(7.51) 


with 
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(7.52) 


(C33bl  +  C13R1) 


'33  2  °13I'2J 


o'  33  2  ''13"2I 

2y1C13 _ 

xoClllC33b2  +  C13R2) 


(7.53) 


(7.54) 


SuDstitution  into  (7.44)  and  (7.45)  yields  the  following  coupled 
equations  for  A*^n  and  C*^: 

TQ4  */» 

Clra^5m  =  "  (1*R2b2}  tanh  tb2jzoi  “J  f  T(x)  cosj*dx 


«  xa 


.  b5  f  aT 

7  Jo  “ 


sinjxdx 


(7.55) 


+  2_  (-l)r 
xo  J 


a2(a2-S-'-  Q*T(xq)  V'  ~r^; 


,2  (-l)m  Y'  /  2  ral{arSl}  Q*  a2(a2‘s2) 

x - 1 -  /  i  I  1  k  — 2~2 - ~ 7? — 2“ 2~" 


xo  J 


t  Z  *  _  Q*  *• .  >  a* 

1  Y7§F\1d 


Aln^4n  ■  -  Q|T(xQ)(a2-S2)  tanh  (a2kxQ) 


„  2b2<1-R2b2)  (-1)"  V 

z  k  Z-f 


k  m  1 


k  +b2j‘ 


*oi  y 

Jq 


aT  .  .  . 

—  sin  jxd* 


JOJ  T(x)  cosjxdx 


(7.56) 


2  (-l)n  V'  (  1xn  -2  F  bi^1“Rfbi>  ,  b2(1-R|b2)] 

77  k  Lj  M)  j  T7T272  ^3  T7TZ7Z  c*n 

o  m  k  +d^j  k  +D.,J 


k  +b2j 
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where 


^5m  =  tann  <W  '  Q$f1-R2b2^  tanh  ^b2^zo^ 

(7.57) 

^4n  =  ^ai-Si^  tanh  (ai^0)  ~  Q*( a2-S2)  tanh  (a2kxQ)  • 

Equations  (7.55)  and  (7.56)  are  also  of  the  form  of  (7.34)  and  can 
be  solved  with  the  same  iteration  procedure. 

7.3  COMPUTATIONAL  RESULTS 

In  this  section  we  present  calculated  results  for  several 
discontinuous  temperature  distributions.  The  input,  plane-stress 
elastic  constants,  tne  same  as  those  used  for  the  calculations 
reported  in  Reference  9,  are  appropriate  to  a  50  percent  aluminum- 
graphite  composite, 

=  3.187  Msi 

=  1.118  Msi 

55.401  Msi 
=  3.298  Msi 

=  5.033  x  10~5  Msi/°F 

=  3.705  x  10"5  Msi/°F. 

expansion  coefficients  are 

=  15.705  x  10"6  °F-1 

=  0.352  x  10"6  °F-1 

and  tne  length  to  width  ratio  is  10, 

*0  ■  1 

zo  *  10 

Figure  7.1  is  a  plot  of  the  shear  strain  as  a  function  of  z 
along  x  =  0  for  the  odd  temperature  distribution 


°13 

C33 

C44 

Y1 

y3 

Tne  tnermal 

el 

e3 
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200  Terms 
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Figure  7.1.  Shear  strain  along  z  =  0  for  the  odd  temperature  distribution  T, (x). 


T±(x)  = 


1 

'  7  * 


"  x0  1  x  <  0 


(7.58) 


+ 


1 

2 


0  <  x  1  x0 


Tne  solid  line  indicates  the  shear  strain  as  calculated  with  200 
coefficients,  and  the  dots  indicate  selected  values  as  calculated 
with  100  coefficients.  The  numerical  results  give  no  indication  of 
a  singularity  in  the  stress  field;  and,  indeed,  from  Bogy's  work  on 
wedges  (Reference  11),  one  would  not  expect  one. 


Figure  7.2  shows  the  results  of  two  calculations  for  the  shear 
strain  as  a  function  of  z  along  x  =  0.8  xQ  for  tne  odd  and  even 
temperature  distributions. 


xQ  <  x  <  -  0.8  xQ 

0.8  x„  <  x  <  +  0.8  xrt 
o  o 

0.8  x„  <  x  <  x„ 
o  o 


(7.59) 


and 


T3(x)  . 


-  x  <  x  <  -  0.8  x„ 
o  o 


-  0.8  xq  <  x  <  +  0.8  xQ 

0.8  x  <  x  <  x 
o  o 


(7.60) 


The  shape  of  the  curve  corresponding  to  the  even  temperature 
distribution  T^(x)  was  unexpected,  but  so  far  no  algebraic  or 
coding  errors  have  been  found.  In  any  case,  it  has  been  verified 
that  the  longitudinal  force  at  z  *  0, 


<7zz(x,0)dx 


> 


i  on 


T,(x) 
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Figure  7.2.  Shear  strain  along  x  =  0.8xo  for  the  odd  and  even  temperature 
distributions  T0(x)  and  T~(x). 


balances  that  of  the  shear  stress, 

Jxz  *°'8  Vz)dz  * 

A  plot  of  the  longitudinal  strain  e2Z  as  a 
along  x  =  0.8  xQ,  corresponding  to  the  temperature 
T2(xj  and  T^(x),  is  given  as  Figure  7.3. 


function  of  z 
di  stributions 
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along  x  =  0.8xo  for  the  odd  and  even  temperature 


8.  SUMMARY  AND  CONCLUSIONS 

a.  For  G/A1  where  the  matrix  is  in  a  T4  temper  (rapid  quench  to 
-240°F) : 

Axial  stresses  in  the  aluminum  at  room  temperature  are  tensile 
and  lie  in  the  range  30  to  40  ksi . 

Subsequent  reheating  will  lower  the  axial  stress,  and  the 
matrix  will  remain  elastic  up  to  a  temperature  of  approxi¬ 
mately  600°F,  at  which  point  the  axial  stress  will  have  become 
compressive. 

Residual  hoop  stresses  in  the  matrix  are  tensile  and  range 
from  approximately  10  to  20  ksi  for  all  cases  in  which  the 
fioers'  transverse  coefficient  of  thermal  expansion  remains 
less  than  or  equal  to  that  of  the  matrix. 

The  interfacial  radial  stress  is  compressive  and  ranges  from  2 
to  10  ksi  so  long  as  the  fibers'  transverse  coefficient  of 
thermal  expansion  does  not  exceed  that  of  the  matrix.  This 
would  appear  to  be  a  most  desirable  result  in  maintaining  the 
integrity  of  these  materials,  since  they  are  currently  limited 
oy  their  aoility  to  bear  significant  radial  tensile  stresses 
across  this  interface. 

o.  For  G/Al  where  the  matrix  is  in  a  TO  temper  (slow  quench  to 
-240°F),  tne  axial  residual  stresses  are  near  zero  at  room 
temperature. 

c.  For  G/Al  in  the  T4  temper,  a  reduction  of  approximately  5  to 
10  ksi  in  the  room  temperature  axial  residual  stress  could  be 
expected  if  the  composite  was  quenched  to  -320°F. 


133 


For  tne  A/k\  and  SiC/Al  composites,  residual  stress  results 
are  qualitatively  similar  to  those  of  G/Al  of  equivalent 
temper. 

Calculated  axial  stress  strain  behavior  o^  SiC/Al  using  a 
yielded  matrix  material  in  the  TO  temper  condition  is  in 
re  sonable  agreement  with  experimental  data  up  to  200  ksi. 

For  a  three-layer  laminate  of  G/Al ,  minimum  longitudinal 
coefficient  of  thermal  expansion  occurs  at  a  ply  angle  of 
approximately  20  degrees  (i.e.,  20/0/20). 

PRUFC  yields  G/Al  thermoelastic  properties  in  good  agreement 
with  those  reported  by  Hashin  and  Humphries  who  used  a  finite 
element  analysis  (Ref. 2  ). 

The  S-C  'BED  residual  stress  code  yields  results  for  G/Al  which 
are  in  good  agreement  with  the  ANSYS  finite  element  results 
reported  by  Hashin  and  Humphries  (Ref.  2). 
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APPENDIX  A 
ALUMINUM  PROPERTIES 

Thennoelastic  properties  and  yield- strength  models  for  .6061 
and  2024  aluminum  alloys  are  summarized  below. 


A. 1  6061  ALUMINUM 

for  the  6061  aluminum  matrix,  the  elastic  modulus  was  input  as 
a  function  of  temperature,  linear  interpolation  being  used  between 
tne  following  tabular  points  Jn  Table  Al: 


TABLE  Al 

ELASTIC  MODULUS  FOR  6061  ALUMINUM 


T(°F) 

-459 

75 

240 

400 

600 

1100 

The  Poisson's  ratio  was  also 
according  to 

v  =  0.5  -  E/58.823 


E(Msi) 

11.8 

10.0 

9.8 

9.0 

7.0 

0.0 

allowed  to  vary  witn  temperature 

( A.  1) 


where  E  is  the  elastic  modulus  in  Msi  as  computed  from  the  table 
aDove.  The  thermal  expansion  coefficient  for  tne  aluminum  matrix 
was  taken  to  be  temperature  dependent  according  to  the  tabulated 
points  in  Taole  A2: 
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TABLE  A2 

LINEAR  COEFFICIENT  OF  THERMAL  EXPANSION  FOR 
AT  6061  AND  A1  2024 

T(°F)  a( 10~6°F_1 ) 

-  400  0 

-  300  5 

-  200  9 

0  12 

1100  18.5 

o 

Residual  stress  calculations  were  done  with  two  different 
yield  strength  models.  A  model  appropriate  to  6061  A1  in  the 

annealed  state,  or  TO  temper,  is  depicted  in  Figure  Al.  The 

room- temperature  yield  stress  at  0.2  percent  offset  if  8  ksi  in 

accord  witn  Ref.  12  (the  elastic  modulus  at  room  temperature  is 
taicen  as  10  Msi).  A  yield  model  appropriate  to  6061  Al  in  the  T4 
temper  (rapid  quencning)  is  shown  in  Figure  A2.  Here  the  yield 
stress  at  room  temperature  at  the  0.2  percent  offset  point  is  21 
ksi,  as  given  in  Ref.  12.  Tne  shapes  of  the  room-temperature 

stress-strain  curves  for  both  tempers  were  scaled  from  that  for  Al 
201  in  tne  0  temper,  for  whicn  experimental  data  were  available  (see 
p.  25  of  Ref.  1).  The  variation  of  yield  stress  with  temperature  is 
scaled  proportionately  to  that  of  other  aluminums  (see,  for  example. 
Ref.  13). 

A. 2  2024  ALUMINUM 

For  tne  2024  aluminum  matrix,  tne  elastic  modulus  as  a  func¬ 
tion  of  temperature  was  given  by  the  tabulated  points  in  Table  A3: 
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r.  no- -  in/in) 

Mgure  A-T.  y,-  ,  , 

temper S(annea Ud)!  ^  6061  aW 


TABLE  A3 

ELASTIC  MODULUS  FOR  2024  ALUMINUM 


T(°F) 

-  459 
75 
240 
400 
600 
1100 

The  Poisson's  ratio  was 
according  to 

v  =  0.5  -  E/61.765  , 


E(Msi) 

12.39 

10.5 

10.29 

9.45 

7.35 

0 

allowed  to  vary  with  temperature 

(A. 2) 


where  E  is  the  elastic  modulus  in  Msi  from  Table  A3.  The  thermal 
expansion  coefficient  was  taken  to  De  the  same  as  tnat  for  A1  6061 
as  given  in  Table  A2. 

The  yield  stress  model  for  tne  TO  (annealed)  and  T4  tempers 
are  shown  in  Figures  A3  and  A4. 
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APPENDIX  B 
LAMINATE  CODE 


A  description  is  given  here  of  the  code  constructed  at  S-CUBED 
for  calculating  the  thermoelastic  properties  of  a  cross-plied 
laminate  wnose  individual  plies  may  be  anisotropic- transversely 
isotropic.  The  input  properties  may  be  either  the  overall  elastic 
constants  and  thermal  expansion  coefficients  of  each  type  of  ply  (up 
to  three  types  oriented  in  three  different  directions)  or  the 
corresponding  quantities  for  the  fiber  and  matrix  constituents  in 
each  type  of  ply. 

The  required  input  deck  is  described  in  Section  Bl,  a  sample 
output  is  given  in  Section  B2,  and  a  FORTRAN  listing  of  the  code 
constitutes  Section  B3  below. 

B.l  INPUT  DATA  DECK 

Any  consistent  set  of  units  may  be  used  for  the  input 
parameters;  the  corresponding  output  will  be  in  the  same  units. 


Card 

No. 

Format 

Columns 

Input 

Para¬ 

meter 

Description 

1 

110 

1-10 

NPROB 

Number  of  calculations 

2 

20A4 

1-80 

HEADER 

User's  problem  descrip¬ 
tion 

3 

8E10.0 

1-10 

VFRAC(J) 

Volume  fraction  of  first 
type  of  ply 

11-20 

ANCLE (J ) 

Orientation  of  first 
type  of  ply  with  respect 
to  laminate  z-axis 
(degrees)  (See  discus¬ 
sion  below). 
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Card 

No.  Format 

Columns 

Input 

Para¬ 

meter 

Description 

21-30 

TYPE(J) 

Type  of  input  to  follow. 

r  jt  t unr  a  •  . .  r-  ^  i_ 

If  TYPE>0  input  for  the 
individual  fiber  and 
matrix  components  is 
expected  (2  additional 
cards) 


o 

If  TYPE<=0,  input  of 
overall  ply  properties 
is  required  (1  addi¬ 
tional  card). 

3A 

8E10.0 

Ply  proper¬ 
ties 

Card  3A  is  required  if 
ply  properties  VFRAC>0 
and  TYPE<=0.  Otherwise 
omit. 

1-10 

ET(J) 

Transverse  elastic  modu¬ 
lus  for  first  type  of  ply 

11-20 

EA(J ) 

Axial  elastic  modulus. 

21-30 

NUTT(J) 

Poisson's  ratio:  trans¬ 
verse  contraction, 
transverse  stress. 

31-40 

NUTA(J) 

Poisson's  ratio:  trans¬ 
verse  contraction,  axial 
stress. 

41-60 

GA(J ) 

Inplane  shear  modulus 

51-60 

ALFT(J) 

Transverse  linear  coef¬ 
ficient  of  thermal  ex- 
pansi on 

61-70 

ALFA(J) 

Axial  linear  coefficient 
of  thermal  expansion. 

3B 

8E10.0 

Fiber 

Properties 

Cards  3B  and  3C  are  re¬ 
quired  if  VFRAC>0  and 
TYPE>0.  Otherwise  omit. 
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Card 


No. 

Format 

Columns 

meter 

Description 

1-10 

ETF(J) 

Transverse  elastic  modu¬ 
lus  of  fiber  in  first 

ply. 

11-20 

EAF(J) 

Axial  elastic  modulus  of 
fiber. 

21-30 

NUTTF(J) 

Trans verse- transverse 
Poisson's  ratio. 

31-40 

NUTAF(J) 

Poisson's  ratio:  trans¬ 
verse  contraction,  axial 
stress. 

41-50 

GAF(J) 

Transverse-axial  shear 

modulus  for  fiber. 

51-60 

ALFTF(J) 

Transverse  linear  coef¬ 
ficient  of  thermal 
expansion. 

61-70 

■ALFAF(J) 

Axial  linear  coefficient 
of  thermal  expansion. 

71-80 

VF(  J ) 

Fiber  volume  fraction  in 
first  type  of  ply. 

3C 

8E10.0 

Matrix 

Identical  to  Card  3B, 
except  that  the  proper¬ 
ties  correspond  to  the 
matrix  material  in  ply  J. 

The  sets  of 

data  cards  2 

,  3,  ...  3C 

are  repeated  as  many  times 

as  specified  on 

Card  1.  For 

each  set  of  data,  Card  3  must  appear 

three 

times,  once  for  each 

of  tnree 

different  orientations  of 

plies. 

However, 

if  VFrtAC  is  set  to  zero 

on  Card  3,  no  further  data 

for  that  type  of  ply  is  required. 


If  the  orientation  ANGLE(J)  is  set  to  0.0  on  all  three  of 
Cards  3,  the  code  will  run  a  series  of  calculations  for  orientations 


+e,  -e,  0.0  degrees  of  the  three  ply  types,  for  increments  in  e  of  5 
degrees  from  0.0  degrees  to  90  degrees.  Otherwise  the  code  will  do 
one  calculation  of  laminate  properties  for  the  ply  orientations 
speci  fied. 

B.2  SAMPLE  OUTPUT 

A  sample  output  is  given  below  for  the  case  of  a  laminate  of 
SiC/Al  plies,  and  in  this  case  properties  of  the  fiber  and  matrix 
constituents  were  input  for  each  type  of  ply. 

The  first  section  of  output  is  a  printing  of  the  user-supplied 
neader  followed  by  a  reiteration  of  the  input  properties  (fibe*  and 
matrix  for  each  ply  type  in  this  case).  The  properties  as 
calculated  with  tne  PRUFC  code  (ply  coordinate  system)  are  then 
given  for  the  three  ply  types,  and  the  stiffness  coefficients  and 
stress-temperature  coupling  coefficients  are  also  printed. 

The  main  output  is  then  begun.  In  tnis  example,  the 
orientation  angles  ANGLE(J),  J=l,2,3,  were  set  to  zero,  so  a  series 
of  output  results  at  +e,  -e,  0.0  degrees  are  printed,  only  one  of 
whicn  is  included  here  for  the  +35,  -35,  0.0  degree  case 
corresponding  to  input  volume  fractions  e  0.25,  0.25,  a. .3  0.50 
respectively.  The  stiffness  and  complice  ~M  . "ficients  are  listed 
followed  Dy  the  thermal  expansion  cjer  i-  ts.  Finally  the 
intraply  stress-tanperature  derivatives  (K-.  r  v*  t  in  this  case)  are 
given  in  Doth  the  laminate  and  ply  coordinate  .yb-ans. 

Tne  final  bit  of  output  is  a  summary  for  each  of  the 
orientations  (5  degree  increments)  of  the  engineering  elastic 
constants  and  tnennal  expansion  coefficients  for  the  laminate. 
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B.3  CODE  LISTING 


A  FORTRAN  listing  of  the  code  is  given  below. 
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APPENDIX  B3 

LAMINATE  MAIN  PROGRAM 

REAL  Nx,N2,N3,NUXZ,NUYZ,NUTT,NUTA,NUTTF,NUTaF,nUTTM, YOTA" 
DIMENSION  S(b,b),SINV(b,6),AX(l9>,AV(l9>.AZ( 19 1,6 AA 11 91 
DIMENSION  HEAOERIZO},  A|*l,  SI600T 0,31,  PHI (31,  SIS»LY(9,3) 
DIMENSION  ClN(8,3),CPRlH£(b,6> .Ell  19) ,E2 1 19 | ,£3 1 1 » > , RUIZ <19 1 , 

.  NUYZ119),6XY119),GXZU9),TYPE(3>»C0UH(8),GPRIMd9) 

C0MM0N/0NE/N1,  N2,  N3 

C0MN0N/TW0/C1)9,9),  C2(9,9),  C3 19,91 
COMMON/THREE/GAMMAK*),  GAMMAJ|9|,  GAMMA 3 19  1 
C0MM0N/F0UR/01(2,2> ,  02(2,21,  03(2,21 

COMMON/INPUT /9FR AC (31 , ANGLE (31 ,eT(3) ,E • ( 3 > , rUT T 1 3),NJTA( 31 ,GAI 
.ALFTt3l,ALFA(3l,ETF(3),CAF(3l,NUTTF(3>,NUTAF(3>,6AF(3l,ALFTF(2 
.ALFAF)3)«vF(3),ETH(3),E  AM 1 3 ) , NUT  T  M 1 3 I ,NUT  AM13)»GAN1J)»ALFTN1JI 
•ALfAM(3I ,VM(3) 

DATA  PI/3.1H15927/ 

REA0(S,101>  NPROB 
NNNSO 

i  continue 

READ'S, 102)  HEADER 
WRITE16, 10DI 
HRITE( 6 , 112)  HEADER 
DO  11  J=l,3 

RE  ADIS, 103)  VFrACIJ),  ANGLE(J),  TYPEfJ) 

IFI¥FRiCIJ).LE»0.»  SO  TO  18 
1F1TYPC) J) .ST.O.)  GO  TO  12 


READ  IS  ,  1 
WRITE ( 6 , 
WRITe(6, 
WRITE) 6* 
<50  TO  11 

18  DO  19  X  = 

19  CINIk.J) 
60  TO  11 

12  R£A0(S,1 
.J),VF< J) 

RE  AO ( S , 1 
.J) ,WH( J) 
WRITE (6, 
WRITdb, 
WRITE) b, 
WRiTOb, 
WRI TE ( b • 
WR I TE ( 6 , 
WRITE(b, 

ii  continue 
NUVFRAC 
N2-VFRAC 
N3=VFRAC 
00  13  JS 
IF ( YFRAC 
IF ) TYPE ( 
CALL  PRO 

13  CONTINUE 

DO  19  JZ 

IF) yFRAC 
IF(TYPE) 
CALL  CCO 

19  CONTINUE 
DO  22  J= 
WRITE) b, 
WRiTOb, 

22  WRITdb, 
WRI TE ) 6 , 
WRiTOb, 
WRITdb, 
ANGL=-5. 
9  MA  X  - 19 
IF ) ) ABS) 
00  10  MK 
AnGL=AMG 
IF ) NMAX ■ 
THETASAN 


W0/C1 19 ,9  )  ,  C2 ( 9 ,9  I  ,  C3I9 ,9 ) 

HrEE/GAMMAI )9 ) ,  GAMMAS  )9  )  ,  GAMMA 3 ) 9  ) 

OUR/01 (2,2),  02(2,2),  03(2,2) 

NPUT/YFRAC(3), ANGLE'S), eT(3),EA(3>,nUTT(3),M3TA(3),GA(3), 
ALFA)3),ETF(3),CAF)3),NUTTF(3),NUTAF(3),6aF)3),ALfTF(3), 
,yF)3),ETM)3),EAM<3),NUTTM)3),NUTAM)3),SAM) j),ALFTM)3), 
,YM(3> 

3.1915927/ 

01)  NPROB 


03)  ET(J),EA)J) ,NUTT( J),NUTA) J) ,GA ) J ) , ALF T ) J) , ALF A ) J I 
10S)  Jr  VFRAC(J),  ANGLEIJ) 

1061  CT)J),  EA)J),  GA)J| 

107)  NUTT(J),  NUTA1J),  ALFT1J),  AlFAIJ) 


l03)  ETF) J),EAF( J),NUTTF) J),NUTAF(J),GAFI J),AlFTF)J),ALFAF) 
) 

103)  ETM(J),EAM( J),NUTTM|J),NUTAMIJ),GAM|J|,ALFTM)J),ALFAMI 
) 

,108)  J,  YFRAC  I J  )  ,  AM6LE '  *)1 
,109)  YF(J) 

,106)  ETF 1 J ) ,  EAFI J) ,  BAFIJ) 

,107)  NUTTF(J),  NUTAF1J),  ALFTFIJ),  ALFAF(J) 

•110)  VMIJ) 

, 10b)  ETN(J),  EAMIJ),  GAN(J) 

,107)  NUTTN(J),  NUTAM(J),  ALFTM1J),  ALFAMIJ) 

E 

cm 

C(2) 

03) 

=  1,3 

C(J).LE.Q.)  GO  TO  13 
(J).LE.O.)  60  TO  13 
UFC) J,CIN) 

E 

=  1,3 

C)J).LE.O.)  GO  TO  19 
(J).GT.O.)  GO  TO  19 
OCF! J,CIN) 

E 

=  1.3 
•  115)  J 

,11b)  )CXN)1,J),I=1,9) 

,117)  )ClN!I,J),I=$,S> 

,100) 

,112)  HEADER 


IF  ( t ABS  I  ANGLE  1 1 ) ) • A  BS I  A  N6LE I  2 ) )♦ ABS) ANGLE) 3)l).GT.O.)  XMAxri 
00  10  NK;1,NMAX 
AnGL=AMGL  *  5. 

IF1XMAX.C0.1I  ANGLIANGLE ID 
TMCT  a: ANGL*P 1/180. 

00  i5  x : l , 8 

IS  COUMfX  CCINIK,  I ) 

CALL  CR0T1CPRIMC .GPRIME ,CDUM , T HE T A ) 

00  20  JSl ,9 
00  20  I : I , 9 
20  C11I,J)=CPRIME1I,J) 

00  21  J:I,9 
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S3* 

21 

GAHHAlf  JirGPBZKE  I  J> 

81* 

00  23  J=l,2 

•  2* 

00  23  1=1,2 

83* 

23 

Oil l,J)=CP»IM£f 1*8, J*8I 

88* 

THETA:  -  THET* 

85* 

IFIKHtX.EO.il  ThETA=AN6lE<2)*PI/18C. 

86* 

00  16  K=1 , 8 

87* 

lb 

C0UH«K)=CIN|K,2> 

88* 

CALL  C(T0TICFI»1HE,GP»1HE  ,CDU1, THETA) 

89* 

00  30  J=1 ,8 

9(j* 

00  30  1=1,8 

91  * 

30 

C2II,J):CF8IHEI1 ,J) 

92* 

00  31  J=l,8 

93* 

31 

GAHHA2IJ 1=6991 HE IJ) 

9** 

00  28  j= 1 , 2 

95* 

00  28  1=1,2 

95* 

28 

02II,JI=CP8IHEII*8, J«8I 

97* 

THETASO. 

98* 

IMKHAx.EO.l  1  THE  TA :ANGLE  |3|*9I/180. 

99* 

00  17  k-1,8 

100* 

17 

C0UHIK>=cINIK,3) 

101* 

CALL  CftOTiePTXHE.GPRlMC ,COUH, THeTA) 

102* 

00  80  J=1 , 8 

103* 

00  80  1=1,8 

108* 

80 

C3II,JI=CpRIHEH,JI 

105* 

00  81  J=1 , 8 

106* 

81 

6AHHA3IJI=GP91HEIJI 

107* 

DO  25  J=l(2 

108* 

00  25  1=1,2 

109* 

25 

03II,JI=Cp8lHEU*8,J*9l 

no* 

CALL  EXIS11,S21,S31,S81I 

111* 

CALL  ETISi2,S22.S32.S<i2) 

112* 

CALL  E2(S13,SZ3fS33*S83) 

113* 

CALL  GTZf S18,S28,S38,S88I 

118* 

CA(.l  GXZXVIsSS, 556. 565,5661 

115* 

00  26  J: 1 , 6 

116* 

00  26  1=1,6 

117* 

26 

S|I,J)=0. 

118* 

SII.IISSU 

119* 

511,21=512 

120* 

511,31=513 

121* 

511 ,8)=2.*S18 

122* 

3I2,1|=S21 

123* 

512,21=522 

128* 

512,31=523 

125* 

512, *1=2. *528 

126* 

SI3,l)=S3i 

127* 

5I3,2)=S32 

128* 

SI3,3I=S33 

129* 

51 3, *1=2. *538 

130* 

SI8, 11=581 

131* 

518,21=582 

132* 

518,31:583 

133* 

SI8,8|:2.aS88 
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13** 

135* 

136* 

137* 

138* 

139* 

1*0* 

1*1* 

1*2* 

1*3* 

x*** 

1*5* 

1*6* 

1*7* 

1*8* 

1*9* 

ISO* 

1S1* 

152* 

153* 

IS** 

155* 

156* 

157* 

158* 

159* 

160* 

161* 

162* 

163* 

16** 

165* 

166* 

167* 

168* 

16»* 

170* 

171* 

172* 

173* 

17** 

175* 

176* 

177* 

178* 

179* 

180* 

181* 

182* 

183* 

18** 

185* 

186* 

187* 


S 1 5, 5) =555*2 . 

S!S, 61=556*2. 

SI  6, 5  I :S65*2  • 

S I  6  «  6 1  =  S66*2 • 

CALL  INVRSI6,S,SInV» 

CALL  THERNLI A.SIGOOT) 

AX  I KK )  =  A ! 1 ) 

Ay I NN I =A I  2  I 
AZIKK)=A!3) 

EltKKl=l./Sll 
C2!KK)=1./S22 
E3IK*I=1./SJ3 
NUXZIKK):  -  S 1 3*E  3 1  NX  t 
NUTZ<KK>=  -  S23*E3IXK) 

GaAIKKI=.5*SINVI*,*I 

6X2INK)=.5*SINVIS,S) 

GXT!KK)=,S*S1NVI6,6> 

AN6L2=  “  ANGL 
ANGl3=0. 

IflKNAX.NE.il  60  TO  51 
AN6L2- ANGLE  I  2 1 
AN6L3:AN6LEl 31 
51  CONTINUE 

WRITCI6.1U)  ANGL,  AN6L2,  ANGL3 
WR1TCI6.1181 

WR I Tf  f  6 , 119  I  SINVII,ll,SlNVtl,2),SINVil,3>,SINVt*,] » 

WRITE  1 6, 121)  SINT! 2,2), SINT! 2,3) ,SINV!«,2> 

WRITCI6.122)  5INVI3,3),SINVI«,3) 

XXS .5*5  ZNV ( « ,*  ) 

XY=.S*S1nV<5,5> 

XTY=.5*STN»IS,6I 

XXY=.S*5INVI(,,6> 

WRITCI6.123)  XX 
WRITEI6.137)  XT,  XTT,  XXT 
WRj  TE 1 6 , 12*  ) 

WRITCt6,125)  511,512,513, SI* 

WRITE  1 6 , 126  >  522,523,52* 

WRITE! 6, 1271  533,53* 

WRl TEI 6,128)  S«* 

WRITE  16, 138)  555 ,S56,S66 
WRITE! 6, 129) 

WRITE  1 6 , 13 1 )  A 
WRI TE 1 6, 135 1 
Oo  *2  J=t,3 

IMSIGOoTI*,J).NE.O.)  GO  TO  *3 
GO  TO  *2 
*3  00  **  1=1,3 

**  If  IA8SISIG00TIZ,  JI/STGOOT!*,  JJ  l.LT.l.E-07)  5 IGOOT 1 1 , J > =0. 
*2  CONTINUE 

RMllll=ANGL*f>l/180. 

PMIl2»rANGL2*PI/18C. 
pm 1 1 JI=ANGL3*PI/180. 

00  62  J=1 , 3 
5 1 =5  IN I PMX | J  M 
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co=cos«pmhji  > 

Sl2=SI**2 
C02=l.  -  SI  2 
SX6PLV<l*J)=SIG00Tt 1,J> 

SIGPLYI  2, JI=SlGD0T<2tJ) *C02  *  2*«S1S00T l* , J)*SI«CO  ♦  SIS00T«3tJI 
.  *SX2 

SIGPLTI 3.JI=SIGD0TI 2, JI*SI2  2. *516001 1  * , J)*SI*C 3  *  SI300TI3,J» 
.  «C02 

SIGPLYI9, JI=SIGDOTI2,JI*5I*CO  ♦  SI600T I * , J )• I  COX  •  SI2I 
.  -  SIGOOTI 3, J)*SI*CO 

62  CONTINUE 
00  60  J=l.  3 

60  WRITE  16  » 136 1  J*  I  SI 600T  1 1 , J) ,1  =  1 » 6 > 

WRITE  1 6 • 161 1 

Oo  61  J=1.3 

61  WRITE ( 6f 136 )  J*  ( SI6PLT  I I*J)*I-l*6l 
WRXTCI6.132I 

10  CONTINUC 

WRITC<6,100> 

WRITCI 6* 1 12)  HEADER 

WRxTC(6*133l 

WRITc'6«13*» 

AN6L=-S* 

00  90  I=1,KNAX 
AN6L;ANGL  ♦  5. 

6N6L2-  -  6N6L 
•N6L3=0. 

IFIKHftX.NC.l)  50  TO  92 
AN6L-ANGLE I 1 ) 

AN6L2-ANGLC I  2 1 
AN6L3=AN6LC| 3) 

92  CONTINUC 

90  WRITE! 6, 130)  AN6L , AN6L2 . *NGC3 ,El 1 1 ) ,C2 1 1 ) , C3 1 X ) ,NUX 21 1)  , NUY2 1 1 ) , 
.  6AA«I),AXII)tAT<I)tA2II) 

WRITE! 6 ( 139) 

AN6L--5 . 

DO  2T  I=1,XHAX 
ANGL -* NGL  ♦  5. 

ANGL2-  -  ANGL 
AN6L3rO. 

IFIKNAx.NE.l )  GO  TO  S3 
ANGl:An6L£ ( 1 ) 

ANGL2=»NGLEI2I 
»NGL 3-ANGLE  I  3 ) 

S3  CONTINUC 

2T  WRX  TE  I  6  »  11  0 1  ANGL,ANGL2tANGL3,GX2II),GXYIII 

nnnsnmn  «  l 

IF INNN.LT.NPROB)  GO  TO  1 
STOP 

100  FORHAT I 1H £ ) 

101  F0PMATII1Q) 

1 32  F qRn A  T  A  2 0 A  N I 
103  F0PN«T(8e10.0) 

ION  F  OPN A  T  t 1 M  *20A« • //  ) 
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*2 

*3 

99 

95 

96 
•  7 

M 

99 

50 

51 

52 

53 
59 

55 

56 

57 

58 

59 

60 
61 
62 
63 
69 

65 

66 

67 

68 

69 

70 

71 

72 

73 
79 

75 

76 

77 

78 

79 

80 
81 
82 
83 
89 

85 

86 

87 

88 

89 

90 

91 

92 
t93 


105  F ORMAT I ///  ,  •  INPUT  PROPERTIES  07  PLY*, 12.*  ¥01  7 7 9 Z- * • 7 7 .5, * 

.  AT*,F7.2,*  DEG  TO  LAMINATE  2  *  **15%/! 

106  FORMAT  t •  ELASTIC  HOoULUS  T9 1 *S VERSE  =  • , Z 1 1 . 5 , ♦  • X I ALs* ,E1 1 .5 , * 

.  AXIAL  SHEAR  MODULUS  =  • ,E 1 1 . 5  I 

107  7 ORMAT f *  POISSONS  RATIO  NUTT • #7jj .8 , •  NUTA  = • . 71 1 . 8 * iOX . • ThE»" 

.  EXP  TRANSVERSE:*, Ell. 5, SX,*AXIAL  =  *»EU. 5) 

108  FORMAT!///, •  CONSTITUENT  INPUT  PROPERTIES  07  PLY  * ,1 2 ,SX , »VOL  7RAC  = 
.*.77.5,»  AT*, FT. 2,*  0E6  TO  LAMINATE  2-AXIS*) 

109  FORMAT  I / »  *  FIBER  PROPER  TIES • ,5* «* VOL  FRAC  =  * ,F7 .5 , / ) 

110  FORMAT!/, *  MATRIX  PROPS RT IES • ,SX ,• VOL  FR AC:*  ,F7.S , / ) 

111  FORMAT!*  PLIES  AT  ♦*,F3.0,*,*,F9.0,*,*,F3.0,*  DEG  TO  LAMINATE  2-ax 
,IS*,/I 

112  FORMAT ( 1M  ,2QA9l 

113  FORMAT!/, 1H  ,9E16.8I 

119  FqRMATIIH  ,8E 19.5) 

115  FORMAT!/, *  STIFFNESS  COEFFICIENTS  OF  PLT*,I2) 

116  FORMAT! IH  ,SX,*Cll=*,£ll.S,SX,*C12:*,Cll.S,5X,*C13=*,E11.5,Sx,’C33 
.:*,EU.S) 

117  FORMAT! 1H  , SX, *C99:* ,£1 1. 5 ,SX , *C66:* »E 11 .5 ,SX , *6A 99 A  1 :* ,E1 1 .S, SX , * 
•GAMMAS:* ,E 11 ,5 1 

118  FORMAT!*  STIFFNESS  COEFFICIENTS  OF  LAMlNATE*l 

119  FORMAT! IH  , SX, • C 1 1  =  * ,El 1 .5 , 5*» *C1 2  =  * ,E1 1 .5 ,5X , *Cl 3= * ,E 1 1  .5 ,5X , *C 19 
,:*,EU.S> 

121  FORhATIIH  ,2SX,*C22:*,£11.S,SX,*C23:*,E11.5,5X,*C29:*,C11.5> 

122  FORHATIIH  ,9SX,*e3S=*,E11.5,5X,*C39:*,ElI.S> 

123  FORHATIIH  ,6SX,*C99=*,Ell,St 

129  FORMAT!/,*  COMPLIANCE  COEFFICIENTS  FOR  LAMINATE*) 

125  FORHATIIH  ,SX,*Sll=*,£ll«5,SX,*S12=*,E11.5,SX,*S13:*,Ell.5,SX,*S19 
.=*,£11,5) 

126  FORHATIIH  , 25X , *522= * ,E 1 1. 5, 5X , *S23= * , El 1 . 5 , 5X , *S29 s • , E 1 1 .5 ! 

127  FORHATIIH  ,9SX,*S33  =  *,En.5,5X,*S39=*,E11.5) 

128  FORHATIIH  ,6SX,*S99=* ,E11.S) 

129  forhati/,*  thermal  expansion  coefficients  for  laminate*) 

131  FORHATIIH  ,SX,*ALF1=*,£11.5,5X,*ALF2=*,E11.5,SX,*ALF3=*,E11.S,5X,* 
•*LF9=»,E11.5I 

120  FORHATIIH  .9E16.8) 

130  FoRMATI*  ♦*,F3.0,*,*,F9.0,*,*,F3.0t9El2.5) 

132  FORMATI/,*  ****•*••****•**««•«••••*••****•••»**••«*•••»•••••*'*,/) 

133  FORMATI//,*  LAMINATE  ENGINEERING  ELASTIC  CONSTANTS  A<«0  COEFFS  OF  T 
.hermal  EXPANSION*) 

139  FORMATI/,*  PLY  ANGLES  * ,5X , *E X • , 1QX , *EY • , 1 pX , *E2* , 9< , • NUX2* , 8X, * 
•NUY2* ,9X,*GY2»,8X,»ALFX»,BX,*ALFT*,8X,*ALF?*,/I 

1J5  FORMATI/,*  PLY  STRESS-TEMPERATURE  DERIVATIVES  •  L,A  MI  NA  TE  COORD  SYST 
•EM)*,/) 

136  FORMATI IH  ,SX,*PLV*,I2,SX,*DSI6XX0T:*,E11.S,SX,*0SI6TT3T=*,E11.5,5 
*X,*0SIGZ20T:*,Ell.S,5x, *OSI GYZDTz* ,E 11 .5 > 

137  FORMiTIlH  ,25X, *C55=* ,C 1 1 • 5 ,5X , *CS6= * ,E 1 1 . 5 ,5X , *C66 : * , E 1 1 . 5 ) 

138  FORMAT! IH  ,25X,*S55=* ,E 11 .5,5X , *S56= * , El 1 .5 , 5X , »S66=* , E 1 1 . 5 ) 

190  FORMAT!*  ♦*,F3.0,*,*,F9.0,*,*,73.0,2E12.5) 

139  FORMAT!/,*  PLY  ANGLES  * , 5 X , *GX 7 • ,9* , *GXT * , / I 

191  FORMATI/,*  PLY  STRCSS-TEMPE RA TURE  DERIVATIVES  «PLY  COORD  SYSTEM)*, 

./! 


9 


ENO 
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1*  SUBROUTINE  CCOEFIJ.CINI 

2*  REAL  NUTTtNUTA,NUATfNUTTF,NUTAFtNUTTM,NUT»M 

3*  DIMENSION  CINIB.3) 

*•  COMMON/ INPUT /VFR AC  I  3|  .ANGLE  (31  ,ET  «  5 »  ,EAI  31  ,  NUTTI  3»  ,  NJMI  31  ,5AI  31  , 

S«  .ALFTI3I*  A|.FAI3>fCTFI3>.EAFI3lfNUTTFI3>  ,NUT»F f J I , 6»F I  Jl  ,»IFTF  (  3  1  . 

A*  .ALFAFI3) «VF(3ltETHI3ltEANt3)»NUTTNl3lfNUrAMI31*6AMT3l«ALFTM(3|* 

7*  .ALFAMI 31 ,¥MI 31 

B*  NUAT=ETI J)*NUTA t J)/EA(J| 

9*  OELTAS1. 

0*  CINll.JI 

I*  CINI2.JI 

2*  C1NI 3f J) 

3*  ClNI A* Jl 

t*  CINIS.JI 

5*  CINI6.JI 

6*  CINC  7, J ) 

7*  CINI8.J1 

a*  return 

9*  END 


END  OF  COMPILATION:  NO  DIAGNOSTICS. 


•  NUTT(JI««2  -  2.PNUTAI J ) *NUA T  -  2. PNUTT I J I *NUT A  I J » PNUAT 
:ETIJI«(I.  -  NUTA( J|*NUAT)/OELTA 
=£Tt JI*|NUTTI Jl  ♦  NUT«IJI*NUAT)/OELTA 
=ET|J)*I1.  ♦  NUTTCJI IANUTAI JI/OELTA 
=EAIJ>*tl.  -  NUTTI J)**2|/DELTA 
=6*« J» 

-«S*(CIN(1«J)  -  CINI2.JM 

-ICINI1.J)  ♦  CINI2.JI IMALFTI Jl  ♦  CINI3«J)*AL?4lj> 
=2.*CIN«3.J|«ALFT|Jt  *  CIN(«,J)*ALF« I Jl 
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1* 
2  * 
3* 
9* 
5* 
4* 
7* 
8* 
9* 
10* 
II* 
12* 
13* 
19* 
15* 
16* 
17* 
18* 
19* 
20* 
21* 
22* 
23* 

29* 

25* 

26* 

27* 

28* 

29* 


C 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  CROTICPRIM£,GPRHE,CIN,THCTAI 
SEE  NCMO  13  OCT  82  RICE  To  SURTMAN 

COMPUTES  CPR INC  1 1 , J I  ELASTIC  CONSTANTS  I*  RoTATEO  C00R3  SYSTEM 

PROM  T ME  CIN  CLASTIC  CONSTANTS  IN  THE  FIBER  SYSTEM 

CINI1)  =  CU,  CINI 2  > :C12 ,  CINI 3 1  ,C1  3  ,  CIN(*l:C3j,  ?INI$)=C99, 

CINI6)=C66  =  .S*IC11  -  C12»  IP  MATERIAL  IS  TRANSVERSE!  Y  ISOTROPIC 

COORO  SYSTEM  ROTaTCO  BY  AN6LE  ThETA  ABOUT  TME  Xl=:ll  AXIS,  WHERE 

THE  FIBER  DIRECTION  IS  THE  OLD  21=31  AXIS 

CINI7)  Z  SAHAX,  Cl N 1 8  I  =  SAMA2 

DIMENSION  CpRlMC 16,61 •  CINI8I,  GPRINC|9» 

X=S INI  THETA) 

y=cos«thctai 

X2=X**2 

Y2rY**2 

X9SX2*X2 

Y9=Y2*Y2 


X3=X2*X 

Y3=Y2*Y 

CPRIMEI 1,1I=CIN<1> 

CPRIMEI 1,2I=Y2*CIN(2I  *  X2*CINI3t 
CPR1MCI1,3I=X2*C1NI2>  •»  Y2*ClNI  3  ) 
CPRIMEI 1,9|=X«Y«<CINI 31  -  CINI 2  I  ) 
CPRIMEI 1,57=0* 


CPRIMEI 1,61=0. 

CPRIM£l2,l»=CPRIMEI 1,21 
CPRIM£I2,2I=YN*C1NI 11  •  X9*CINI9>  * 
.  «  9.*X2*Y2*CIN|5 I 

CpffIMCI2,3l=x2*Y2*IClNll)  ♦  CIN | 9 ) ) 
.  -  9.*X2*Y2*CINISJ 


2.*X?*Y2*CINt  31 
«  1X9  •  Y9l*CINIj| 


30* 

31* 

32* 

33* 

39* 

35* 

36* 

37* 

38* 

39* 

90* 

91* 

92* 

93* 

99* 

95* 

96* 

97* 

98* 

99* 

50* 

SI* 

52* 

S3* 

59* 

55* 

56* 

57* 

58* 

59* 

60* 

61* 

62* 

63* 

69* 

65* 

66* 


CPRIMEI2,9I=  -  X *V3*CIN III  ♦  I X*Y3  -  X3*Y»*CINI3)  ♦  X3*Y*CINI9> 
.  ♦  2«*( X*V3  -  X3* Y | *C IN  I  5  I 

CPRIMC I  2 , S I =0. 

CPR1ME 12,61=0. 

CPRIMEI3,ll=CPRIMCIl,3» 

CPRIMC | 3 ,2 l=CPRIMCI 2,3) 

C PRIME I  3 , 3  I = Y9 #C IN  I  9  |  «  X9*CIN 111  ♦  2»**2*Y2*CINI 3 1  ♦  9.*X2*Y2 
.  *CINIS) 

CPRIMEI 3,9 1=  -  X  3*Y*CIN  111  •  X*Y3*CINI*»  •  IX3*Y  -  X*V3I«CINI3| 
.  •  2.*IX3*Y  -  X*Y 3 1 *C INI  5 1 

CPRIMC I  3,5 1=0. 

CPRIMEI 3,61=0. 

CPRIMCI9,l»=CPRlMCU,9» 

CPRIMEU,2»=CPRIMEI2,9» 

CPRIMCI9,3>=CPRIMCI 3,9) 

CPRIMC |9,9)=X2*Y2*ICINI II  ♦  CINI9II  -  2.*X2*Y2*C I Nl 3 » 

.  ♦  IIY2  -  X2 1 **2 1 *C IN  1 5 1 

CPRIMEI 9,51=0. 

CPRIMCI9«6I=0. 

CPRIMCl5,ll=0. 

CPRIMEI 5,27=0. 

CPRIHCIS,3I=0. 

CPRIMC 15,91 =0. 

CPRIH£tS,5>=Y2«CINl5l  ♦  X2*CINI6> 

CPRIMC |5,6I  =  X*Y*ICINI5I  -  CINI 6 1 1 
C PRIME  16,11  =  0. 

CPR1MEI6,2>=0. 

CPRIMEI6,3)=0. 

CPRIMC 16,91=0. 

CPRIMC|6,5»=CPRIMCI5,6» 

CPRIMC 16, 6I=V2«CINI6|  «  X2*C I N 1 5  I 

SPRiMEU»=CInI7> 

SPRIMCI 2l=V2«cINI 71  ♦  X2«CINI8I 
6PRIMCI3I=X2*CINI7I  ♦  Y  2*C I N I  8 1 
GPRIMCI9|:tCIN|S>  -  CINI7»I*X*Y 
RETURN 
END 


CNO  OF  COMPILATION!  NO  DIAGNOSTICS. 
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1* 

2* 

3* 

9* 

5* 

4* 

7* 

8* 

9* 

10* 

11* 

12* 

13* 

11* 

IS* 

It* 

17* 

IS* 

19* 

20* 

21* 

22* 

23* 

29* 

25* 

26* 

27* 

2S* 

29* 

30* 

31* 

32* 

33* 

39* 

35* 

36« 

37* 

38* 

39* 

90* 

91* 

92* 

93* 

99* 

9S* 

96* 

97* 

98* 

99* 

so* 

Si* 

52* 

S3* 

59* 

55* 

54* 


SUBROUTINE  EMS11,S21,S31,S91) 

DIMENSION  XI 4 ) 

REAL  Nl,  N2.  N3,  K 
COMMON /ONE /Nli  N2 •  N3 
CQMMON/TUO/C 1 1  4* 9) «  C2I 9*91 1  0319,91 
C  OMMON/COE  F/RI4»7) 

NI1,1I=C1I1,1I 
Ml, 21=0. 

Kl Jt3>=0. 

Ml  1,91  =  011 ],2> 

Mil ,SI=C1|1 t3| 

MI1,4I=2.*C1U,9I 
Nl 1,71=1. 

N I  2, 1 1  =  0. 

N I  2,2>=C2< I • 1 1 
MI2,3)=3. 

MI2,*I=C2( 1,21 
Mi?,si=r:< i,3» 

N I  2 , 4 >  =  2  *C21 1,91 
M 12 ,71  =  1 , 

M I  3, 1 1  =  0 . 

Ml  3,21  =  0. 

M I  3, 3I  =  C31 1,1) 

HI 3,9I=C3I 1,21 
Ml 3,SI=C3I 1,31 
MIS,  41:2.  *C3|1, 9.1 
Ml  3,71  =  1. 
kI9,1I=n1*C1I1,2» 

K|9,2I=N2*C2H,2I 
M I  9, 31 =N3*C31 1,2 ) 

M|9,9)=Nl*Cil2«2>  ♦  N2*C2I2,2»  ♦  N3*c3l2.2» 
M|9,5>:N1«C1|2,3I  ♦  N2*C2 12,3)  ♦  NJ*Cjl2,3» 
M|9,61=2.«IN1*CH2,9I  ♦  N2«C2 1 2,9  I  *  M3*C3I2.9II 
M I  9 , 71 =0. 

MIS,il=Nl*ClU»3 1 
M|S,2I  =  N2«C2H,3I 
MIS,3I=N3«C3I1,3| 

nIS,9I=n1*C1I2,3I  ♦  M2*C2 12,31  •  N3*C3|2,3I 
MI5,5l=Nl*cl(3,3l  ♦  N2*C2 1 3,31  ♦  N3*C3I3.3I 
MIS.4I  =  2.*IN1*CH3,9I  •  N2*C21 3,91  •  N3*C3<3,9I> 
Ml S, 7)=0. 

MI4,1I=N1«CU1,9I 
NI4,2)=N2*C2I 1,91 
NU,3I=N3*C3I1,9I 

MI4,9l=Nl*Cll2,«>  ♦  N2*C2I 2,91  ♦  N3*C3 1 2,91 
MI4,SI=N1*C1|3,9|  ♦  M2*C2 13,91  ♦ 

MI4,41  =  2.*IN1*C1I9,9I  *  N2*C2 1 9 ,9 |  .  N3*C3I9,«II 
M I  4, 71  =  0. 

CALL  ROOT  1 4, X ) 

EPSX  =N1 *X ill  ♦  N2*X 12)  ♦  N3*X I j I 

Sll=EPSX 

S21=XI9I 

S31=*<5I 

S9t=XI4l 

RETURN 

END 


£13  07  COMPILATION! 


NO  DIAGNOSTICS. 


I  61 


0 

0 


SUBROUTINE  £V<S12,S22,S32,S12I 

OINENSION  X ( &  I 

REAL  NI(  N2 «  NJ,  K 

COMNON/ONC/NI,  N2,  N3 

COMMON/TwO/C 1 1 4  t A  1 1  C2|«,1I,  C3I1.1I 

C0NH0N/C0EF/M6.7) 


Mi,i):cm,n 

mi. 2i=o. 


M 1 , 3 1 =0 . 

Ml, II  =C  III, 31 


Ml,5l=2.iCl(l,1> 

M 1 ,6l =  -  cm, 21 


M2, 11  =  0. 
K|-,2i  =  C-»|I,ll 


M2, 31=0. 
M2,H=C2< 1,31 
M2, 51=2. lC2fl,1) 


M2, 61  =  -  C 2 1  1,21 
M3, 11  =  0. 


K  13 
K  I  3 
M3 
M3 


21=0. 

3>=C3I 1,11 
ll=C3(l,3l 
5>=2.*C3( 1,41 


21 

25 

30 

iQ 


M3 

Mil 

Ml 

MU 

Nil 

M* 

KM 
M5 
MS 
M5 
M5 
H  I  5 
M5 


61=  -  C3I..21 
1I  =  NHC1<1.3I 
2l=N2«C2l 1,31 
3l=N3*C3l 1,31 
1 1  =  N1*C  j I j, 3  I 


N3*c3 13,31 

N3*C 313,1*1 
N3*C3I2,3I1 


N2*C2 13,31 

Sl=2.i|NliCl(3,ll  ♦  N2*C2I3,*1 
61=  -  (NliCl (2,31  «  N2iC2 12,3) 

n=Ni*cm,«i 

2l=N2*C2« 1,61 

3 1 =N3*C3( 1 ,6  > 

1 I =N liC  1  ( 3 ,1 1  ♦  F2*C2(  3,11  ♦  N31C3I3,«1 
51-  2. *<Nl *0111,11  *  N2*C2 ( 1,11  *  N3*C  3  1 1  ,1 1  I 
61=  -  INl*cll2,«l  *  N21C2I2.1I  •  N3«C3(2.1II 
CALL  ROOT I  5, X 1 

516T 1;C1 1 1 12  l*X  1 1 1  •  Cl  12,21  *  C](2,jl*X(Hl 
SI6V2=C2<1,2I*X(2I  *  C2 ( 2, 2 1  *  C2(2,3HX(1l 
SICV3  =  c3I1, 2HXI3I  •  C3I2.21  •  C3(2,3Hx(1l 
SI6V:niiSI6V1  «  N21SIGT2  *  N3*SjGT3 
EPSX=N1*X(1)  «  N21XI2>  *  N31XI31 
312=EPSX/SIG V 
S22= 1 •/ 5 I6T 


2. *C1(2. 111X151 
2.*C2 1 2, 1 1 ix • 5  * 
2  -  *C3 | 2 , 1 1 *X 1 5  I 


5 

6 
7 

I 


S32=X I i I /S 1ST 
S12  =  XI 5  I /S1GY 
RETURN 
ENO 


EN3  OF  CQNP IL  *  TI ON  i  NO  DIAGNOSTICS. 
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1* 
2* 
3* 
9* 
5* 
6* 
7* 
8* 
9* 
10* 
11* 
12* 
13* 
19* 
15* 
16* 
17* 
18* 
19* 
20* 
21* 
22* 
23* 
29* 
25* 
26* 
27* 
28* 
29* 
30* 
31* 
32* 
33* 
39* 
35* 
38* 
3  *  • 
38* 
39* 
90* 
91* 
92* 
93* 
99* 
95* 
98* 
97* 
98* 


SUBROUTINE  £Z(S13,S23,S33,S931 

DIMENSION  V ( 8  I 

REAL  Nl,  N2,  N3  •  K 

COMMON/ ONE /N 1 ,  N2,  N3 

COMMON /T Wo /C 1(9,91,  C2(9,9l,  C3(9,91 

COMMON /COE 7 /K (6,7) 

9(1,11=0(1,11 
K( 1,21=0. 

K( 1, 31=0. 

M(1,9I=C1( 1,2) 

K ( 1 , 5 1  =  2 .*C1 (1,91 

M(i,8i=  *  cm,]) 

9(2,11=0. 

9(2,2|=C2| 1,11 
9(2,31=0. 
k(2,9I=C2( 1,21 
K(?,5I =2.«C2 (1,91 
9(2,61=  *  C2( 1,31 
9(3,11=0. 

9(3,2|=0. 

9(3,31=03(1,1) 

9 ( 3,9I=C3( 1,21 
9(3, 51=2. *C3(1,9I 
9(3,61=  -  C3 ( 1 , 3  > 

9(9,1)=N1*C1(1,2) 

K ( 9, 21 =n2*C2 ( 1,21 
9(9,3I=N3*C3(1,2> 

9(6, 9l=N]«C 1(2,21  ♦  N2*C2(2.2I  ♦  n3*C3(2.2I 
9(9, 51=2. *(Nl*Cl(2,9»  ♦  N2*C2 ( 2,91  ♦  N3*C3(2,9)1 
9(9,6|=  -  (N1*C1I2,3I  ♦  N2*C2(2,3»  *  N3*C3(2.3M 
M(5.11=N1«C1(1,9I 
9 (5,21=N2*C2< 1,91 
9(5,3l=N3*C3(l,91 

9(5,9|=Nl«Ci(2,9l  ♦  N2*C2(2,91  *  N3*c3 12.9) 
9(5,5ls2.*(Nl=Cl(9,91  ♦  N2*CZ (9,91  ♦  N3*C3(9,9ll 
9(5,61=  -  ( Nl*Cl 13,91  ♦  N2*C2( 3,9  )  ♦  N3*C3(3,9)I 
CALL  ROOT (5,91 

SI6Zl  =  Clll, 31*9(11  ♦  CK2, 31*9(91  ♦  C1(3,3I  ♦  2. 

S16Z2=C2(1 ,31*9(2)  ♦  C2(2,3)*X(6>  ♦  C2I3.31  ♦  2. 

SIGZ3=C3(1,3>*9(31  ♦  C3(2, 31*9(91  ♦  C3(3,31  ♦  2. 

S IGZ=N 1 *Sj sZ 1  ♦  m2*SIGZ2  ♦  N3*SI6Z3 

C  PS  9  =N 1*9 (  1 1  ♦  N2*9 121  ♦  93*9(31 

S13=EPSX/SIGZ 

523=9(91 /51S2 

S33= 1 •/ SIS  7 

S93=9(51/SISZ 

RETURN 

CNO 


END  OF  COMPILATION: 


NO  DIAGNOSTICS. 


*:i 13,91*9151 
*C2«  3, 91*9(51 
•C3( 3.91*9(51 
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1* 

2* 

3* 

»♦ 

5* 

6* 

7* 

8* 

9* 

10« 

n* 

12* 

13* 

18* 

15* 

16* 

17* 

18* 

19* 

20* 

21* 

22* 

23* 

28* 

25* 

26* 

27* 

28* 

29* 

30* 

31* 

32* 

33* 

38* 

35* 

36* 

37* 

38* 

39* 

80* 

81* 

82* 


SUBROUTINE  6X2XY(SS5,S56,S6S,S66> 
DIMENSION  X I  6  ) 

REAL  N1,N2,N3,8 
COMHON/ONE/N1 ,  N2,  N3 

COMNON/FOUR/0 1 1  2  «2 ) t  02(2,21,  D3(2,2> 
C0HH0N/C0EF/8 (6*71 
00  10  J=l«  7 
00  10  1=1.6 
10  8(I,J)=0. 

K(l, 11:01(1.11 

Nl 1,21:011 1.21 

«(1,7>=.S 
K 1 2, 31=021 1,1) 

8(?,8|=0?( 1,21 
8(2,71=. 5 
H(3, 51=03(1, II 
8(3,61=0311,21 
8(3,71=. 5 
K(8,ll=01l I. 2> 

8(8,21=01(2,21 

8(5,31=02(1,21 

8(5,81=02(2,21 

8(6,51=03(1,?) 

8(6,61=03(2,21 
CALL  RQOT ( 6 , X ) 

EPS  1 3=n1*X ( 1 1  «  N2* x ( 3 1  ♦  N3*X(5| 

EPS  I ?=N 1 *X ( 2  I  ♦  N2*X (81  ♦  n3*X (6  I 

S55=,5*EPSl3 

565=.S*EPS12 

8(1,71=0. 

Kl2,7)=0. 

8(3,71=0. 

8(8,71=. 5 
k(S,7I=.S 
8(6,7)=. 5 
CALL  Ro03(6,X) 

EP513  =  N1*X( 1 1  ♦  N2*X I  3 )  ♦  n3*X 1 5  I 

EP$1Z=NI*X( 21  ♦  N2*X ( 8 )  «  N3*l I  6  > 

556= .5«EPS 1 3 

S66=.5*EPS12 

RETURN 

ENO 


E NO  OF  COMPILATION:  NO  DIAGNOSTICS. 
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X* 

2* 

3* 

8* 

5* 

b* 

7* 

M 

9* 

10* 

IX* 

12* 

13* 

18* 

IS* 

16* 

17* 

18* 

19* 

20* 

21* 

22* 

23* 

28* 

25* 

2b* 

27* 

28* 

29* 

3q* 

31* 

32* 

33* 

3b* 

35* 

3b* 

37* 

38* 

39* 

*0* 

*1* 

92* 

*3* 

88* 

85* 

8b* 

*7* 

*8* 


SUBROUTINE  872 1 S 1 8, S28,  S3*  ,  SM  I 
DIMENSION  X(6) 

RE«L  NX,  N2,  N3,  8 

COMMON/ ONE /Nit  N2,  N3 

COMMON/ TWO /C 1  1 8  t8 I #  C2(8,8>,  C3I8.8I 

COMMON/ COE F/K 1 6t  71 

Ktl.llsCill.tl 

8(1,2)=0. 

Mil, 31=0. 

8(1,8>=C1(1,2> 

8(1,S)=C1(1,3> 

1(11,6):  -  2.*cl(1.8) 

8(2,11=0. 

8 I 2,2)=C2( 1 t 1 ) 

8(2,31=0* 

8(2,8)=C2( X,2> 

8(2,St=C2( 1.31 
8(2,61=  -  2.*c2<lf8| 

KI3,lirQ. 

8(3,2)=0. 

8(3,3) =C 3(1,1) 

8( 3,8)=C3( 1,21 
8(3,S>=C3(1,3> 

8(3,61=  -  2«*C3( 1 ,8 l 
M(6,1)=NX«CX(X,2) 

H(8,2):N2*C2(1,2I 
K(8,3)=N3*C3( 1,2) 

8(8,*)=N1*C1(2,2)  ♦  N2*C2 (2,2)  «  n3*C3(2,2) 

8(8,5)=N1«Ci(2,3>  ♦  N2*C2(2,3>  ♦  N3*C3(2,3> 

8(8,61:  -  2.*(N1«C1 (2,8 )  ♦  *2*C2 (2,8)  ♦  N3*C3(2,8>) 
8(S,1)=N1*CX(X,3) 

8(5,2)=N2*C2I1,3) 

8(5, 3)=M3*C3( 1,3) 

8 (5,8)=Nl*Ci(2,3>  *  N2*C2 ( 2 , 3 )  ♦  N3*c3(2.3) 

8(5,S):N1«CX (3,3)  ♦  N2*C2(3,3>  *  N3*C3<3»3> 

8(5,6)=  -  2,*(NX*CX( 3,8)  •  N2*C2 ( 3,8 )  •  N3*C3(3,9I) 

CALL  ROOKS, X) 

SI6YZ1  =  C1(1,8)*X(1)  ♦  C1(2,8)*X18)  ♦  C1(3,8)*X(5>  *  2.*CK8,8) 

SIGV?2=C2(1,8)*X(2>  ♦  C2({,«)*X(8)  ♦  C2(3,8)*x(SI  ♦  2.*C2(8,8) 

S16T23=C3(1,8)«x(3)  ♦  C3(2,8)*X(8)  *  C3(3,8)«X(5)  ♦  2,*C3(9,8) 

Sl6V7=NX*SlG7ZX  ♦  N2*sISTZ2  «  N3*SIGVZ3 
EPSX  =N 1 *X ( 1 )  «  N2«X ( 2 )  ♦  N3*X ( 3  ) 

SX9=.S«EPSX/S16YZ 

S28:.5*X(8  )/SI6YZ 

S39=.S*x<S)/S16vZ 

S88= ,S/S I6YZ 

RETURN 

ENO 


ENO  OF  COMPILATION 


NO  DIAGNOSTICS 


X* 

2* 

3* 

l» 

5* 

6* 

7* 

a* 

9« 

10* 

11« 

12* 

13* 

1** 

15* 

16* 

17* 


C 


SUBROUTINE  INVERT 

INVERTS  3X3  MATRIX 

COMMON /NIX  2/Cl  3,3),  0131,  XX 


COMMON/ SOU T /S I  3, 3 ) 

0=C 1 1, 1 1 «t  C , 2, 2 >*C I  3 , 3)  -  C(3,2l*Cf2,3) l  - 
.  -  CI5-1I«CI2,3II  ♦  Cll,3)*ICI2,l)*C|3,2l 

Sll,l)=lcl2,2)*CI3,3)  -  C I  3,2) *C 1 2 ,3  I  I/O 
StiiZC*  ICI I,2I*CI 3,31  -  CI3,2)«CII, 311/0 
S|1,3I=|C< X,2)*C<2,3>  -  Cl2,?>*CIX,3>)/0 
SI  2, 1 1  =  -  I  Cl  2, 1 1 *C( 3 , 3 1  -  CI3,X>*CI2,3ll/0 
SI2.2I-ICI  1,11  *C  13, 31  -  CI3,1)«CI1,3))/D 
SI2,3)=-  ICI 1,1>«CI2,3>  -  CI2, XI *C 11,311/0 
SI3, j)=ICI 2, 1I*CI 3,2>  -  CI3,1I*CI2,2I I/O 
S( 3,2>=-  ICIX,tl*C(3,2l  •  Cl], XI*CtX, 211/0 
SI3,3I=ICI 1,1I«CI2,2)  -  CI2,1I*CIX, 2)1/0 
RETURN 
ENO 


CIX,2I*ICI2,X)«CI3,3) 
•  CIS, X)*? 12, 2)  I 


IN VOOXOC 
IN V0020G 
IN V0030C 
INvOONOG 

invoosoo 

INVC36O0 
IN  V00700 
INVOOSOO 
INV00900 
IN VQ100Q 

iNvoixoa 

INVQ120Q 
I N  V01 300 
INVOINOO 
INVOXSOO 
INVOlbOO 
IN V0X700 


ENO  OF  COMPILATION:  NO  DIAGNOSTICS. 


1* 

2* 

3* 

M 

5* 

6* 

7* 

a* 

9* 
10* 
XX* 
12* 
13* 
XV* 
15* 
lb* 
X  7* 
18* 


SUBROUTINE  INVRSIN.XlN, XOUT I 
C  COMPUTES  INVERSE  XOUT  OP  INPUT  MATRIX  XIN 

DIMENSION  C I b, 7 ) ,  XIN, 6, 6),  XOUTlb.bl,  Xlb) 
COMMON/COEP/C 
00  30  J21,b 
00  30  1=1,6 
30  Cl  I, JlrxlNiI ,JI 

00  XO  J=1,N 

00  20  K=1,N 
CIK,N«1I=0. 

IFIK.EO.J)  CIK,N*1)=1. 

20  CONTINUE 

CALL  ROOT  I N  t  X  ) 

00  NO  L=i,N 
NO  XOUTIL,J):X|L) 

XO  CONTINUE 
RETURN 
ENO 


END  OP  COMPILATION 


NO  DIAGNOSTICS 


I* 

SUBROUTINE  NIX 

MIXOOIOO 

2* 

C 

COMPUTE  EFFECTIVE  ELASTIC  CONSTANTS  PROM  THOSE  OP  CONSTITUENTS, 

MIX00200 

3* 

C 

Cll  *  C12,  C 13,  C  33,  6AMHA1,  6AMN*3 

MI *00300 

9* 

REAL  Nl,  N2,  K  1 2 1  ,  KE 

Mi X00900 

5* 

CONNON/ELAST/Nl,  Cl 1 1  2  )  ,  C12I2I 

Ml XQ0500 

t« 

COMMON/MIX 1/cll 3,31,  C2I3.3I,  GAMMA  1 (31,  6ANMA2I3),  3ETA1|3|, 

MI *00600 

T* 

.  BETA2I3) 

MI *00700 

8* 

COMMON/ MI X2/CE 13,3),  6AMMAE 13),  KE 

MlXOOSOo 

9* 

Kin -ci (i,i)  •  cm, 2) 

Ml *00900 

io* 

K 1 21 =C2( 1,11  ♦  C21 1 ,2  1 

MIXoiOOO 

u* 

N2- 1 •  -  Nl 

NIX01100 

12* 

H=N2*KI11  «  Ni*K 1 2 1  ♦  C2ll,l)  -  C2I1.2I 

MI *01200 

13* 

KE-N1*M 11)  •  N2«K (21  -  (Kin  -  A ( 2 1 1 •• 2*N 1*N 2/M 

MIX01300 

19* 

CEU,3I=N1*C1(1,3I  ♦  N2*C2 1 1,3)  •  (Kill  -  K  1 2 1 1*1  Cl  1 1 , 3 1  - 

NIX0190C 

15* 

.  C2(1«3))*N1*N2/N 

MI *01 500 

16* 

CE 1  3,31  SNi*Cl  1  3,31  *  N2*C2 1  3,3)  -  2.*IC1U,3)  -  C2I1.31 1**2*N1*» 

12/MI *01600 

17* 

.  H 

HIX01700 

18* 

GAMHAEI 1I=N1*GAHMA1 ( 1 1  ♦  N2*SAMMA2 111  -  1  GAMMA  1 1  1 1  -  5AHMA2I11) 

MI *01800 

19* 

.  •  1 K 1 11  -  K  1  2  I  1  *N1*N2/H 

MIX01900 

20* 

6AMHAEI  :i:Nl*G»MM*i  |J|  «  N2*gAMMa2(3I  -  2,  *(  Cl  (1,3>  *  C2U,3D« 

MI *02000 

21* 

•  I 6AHNA 1 1 1 1  -  GAMMA  2 1 1) 1 *N1*N2/H 

M1X02100 

22* 

CE 1 2,3> :CE 1 1 , 3 1 

MI *02200 

23* 

CEI 3, 1 ) =CE (1,3) 

MIX02300 

28* 

CE ( 3«2>?C£ 11,3) 

MIX02M00 

25* 

68 MM A E 12) 55AHMAC 11) 

MI *02500 

26* 

RETURN 

*i X02600 

27* 

END 

MIX02700 

CNO  OF  COMPILATION! 


NO  01*6N0STICS 


20 

21 


C 


SUBROUTINE  MIXCC 

CALCULATE  COMPOSITE  Cll  ANO  C12  SEPARATELY 
REAL  N1 •  K  (  6 , 7 1  ,  KE 
DIMENSION  X I  6 ) 

COMMON /EL  AST /Nl(  01(21,  C12I2I 

COMMON/ MIX 1 /Cl (3,3),  C2 ( 3 • 3  I •  G ANNA  I ( 3  I  •  GaNN*2(31,  3£T*1(3>, 
.  BET* 2  I  3  I 

C0MM0N/MIX2/CE ( 3  »  3 ) •  GA MMAE ( 3  I (  KE 

common/coef/k 

CllHi:Cll  1,11 
C11I2I:C2I 1,1) 

C12I 1 )sCl | 1 ,2) 

C12t2)-C2( 1,21 
Ip(Nl.LE.OI  Nl* 1 «E“Q5 
CALL  SETUP ( 1 1 
CALL  ROOT  I b,X  I 

SUM:X|1I  «  X  1 2 1  ♦  X  (  3 1  •  XMI 
EtA^C 12(21 /C 11(21 

H  =  -X ( 1 1  -  . 5* ( 3.  «  ETAI*X(2»/£TA  -  .5*11.  -  ETA)*X|3>  ♦  X(M> 


MXC001 
MX  c002 
MXC003 
MXCOOS 
MXCOOS 
MXC006 

mxcoot 

MXCOOS 
MX  COO 9 
MXC010 
MXC011 
MX  CO 12 
MXC013 


MxCOlS 

MXCOlfc 

MXC017 

MXCOlS 


SUM  =  SUM/3.  -  2  ,*H/3 . 

Yiri./SUH 

IF(N1.LE.01  Nl=l.E-05 
CALL  SETUP (21 
CALL  ROOT ( 6 , X ) 

SIGR0:C11(2I*(X(  II  ♦  3.*X(2l  -  X(]l  -  3.*x(HI  -  C12(2I 
C66=  .S«(C11(2I  -  C12(2l, 

SI6RT0=C6fc*(-X(ll  -  !.$*I3.  ♦  ETA l*x < 2 ' /CT*  ♦  ,S«|1.  -  E  T  A  I  *  X  I  3  I 

.  -  3.«(M»  -  1.1 

Y2-S IGRO/3 ,  -  2.*S I GR TO/ 3 • 

Y=.5*(Y1  ♦  Y2 I 
CE( 1,1IS.S*( Y  ♦  KE I 
CE( 1,2I=.5*(KE  -  Y l 
CE(2,2I=CEU,1> 

CE ( 2 , 1 i:CE ( 1,2  • 

RETURN 

ENO 


MXC021 

MXC022 

MXC023 

MXC02M 

MXC025 

MXC026 


ENO  OF  COMPILATION:  NO  0IA6N0STICS. 
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1* 

SUBROUTINE  PRUFCU.CINl 

2* 

REAL  NUXY(2I,  NUXZ (  2  1  t  NUZXI2I,  Hi 

RRU009Q0 

3* 

RE*l.  NUTT.NUT* ,NUTTF  ,NUTAF  ,9UTTN,NUT»N 

M 

REAL  NUC  XYC • NUCXZC  t  NUCZXC 

5* 

REAL  KX(2>,  KZI2I,  MxC.  N  ZE •  ME 

PRUOObOQ 

6* 

OINENSION  Cj N ( 8 , 3 ) 

7* 

DIMENSION  SETAE  1 31 f  EXT2I,  2  > •  Rh012»,  6(21 

•  * 

COMHON/HIXl/Cl(3«3l«  C2  <  3* 3  )  t  6AMMA 1 ( 3  T •  G»NM»Z(3I,  3£T»1<3>, 

RRU01203 

9* 

.  BETAZI3I 

PRU01 300 

10 • 

COMMON/ HIX2/CE 1 3 *3) «  GAMNAEI3I,  ME 

PRU0190Q 

11* 

COMMON/ELAST/Nl,cll(2>t  C12I2I 

PRU01500 

12* 

COMMON /SOUT/S I  3  «  3 1 

PRU01 600 

1J* 

COMMON/ INPUT/ VFRAC (3) .AN6LE ( 3 1 ,E T I  3 1 1 £ A ( 3) . NUT T ( 31 , N JT A I 31 , S»C 31, 

11* 

•AlFTISI ,ALF»t3I.CTF(3l,EAF(3l,NUTTF|3|,NUTAF(3) 

,6AF(3l,ALFTF(3l, 

15* 

.ALFaFI3I,Vf(3I,ETNI3I,EAm(3>*9UTTN(3>,NUTAM(3I, 

SAM( 3| ,ALFTM(  31 , 

16* 

.ALFAMI3I ,VMI 31 

17* 

EX ( 1 |:CTF • J| 

IS* 

EZtll:EAF( Jl 

19* 

NUXV(1I=nOTTF( J) 

20* 

NUXZI1I:nuT«F( ji 

21* 

BETAMl)tALFTF(J) 

22* 

BETA1(3)=ALFAF( Jl 

23* 

KXI1U1. 

29* 

M2 ( 1 1- 1 • 

25*  • 

RHO( l l : 1 • 

2S* 

6(1) -6AF 1 J 1 

27* 

EX(2l=£TM( Jl 

21* 

EZI2I=EAMI J) 

29* 

NUXVI2I=nUTTM(JI 

30* 

NUXZ(2>:NUrxMI Jl 

31* 

BETA2m:ALFTM(J) 

32* 

S£TA2(3|=ALFAH( Jl 

33* 

NX(2I=1. 

39* 

KZ( 2>-l • 

35* 

RM0(2I=1. 

36* 

6(2I=GAH( J) 

37* 

C 

CALCULATE  CSUBI J*S  FOR  FIBER 

PR UO 3 9 00 

31* 

NUZXf li:(Ex(  1  l*NUXZ(l  ii/EZfl  1 

PRU03500 

39* 

0ELTA=1.  -  NUX Y 111 **2  -  2.*NUXZ( 1 )*NUZXI 1)  -  2* 

•  N JX Y 1 1 1 *NU  X  Z (11 

PRU03600 

93* 

.  *NUZX (It 

PRU03700 

91* 

Cl ( 1 • i I :EX I  I 1 • I  1 .  -  NUXZ(1I*NUZXI1) I/DELTA 

PRU03B03 

92* 

Cl(l,2i=EX(l l*(NUXV(l  )  ♦  NUX Z ( 1 )*NUZX (lll/OELTA 

PRU03900 

93* 

CKlt3i:CX(ll«(l.  «  NUXV(1||*NUXZ(1)/0ELTA 

PRUO*  000 

99* 

Cl<3'3)=CZ(l>*(l>  -  NUXY(li«*2l/0ELTA 

PRUO* loo 

95* 

Cl(2»H=CKl,2l 

PRT09200 

96* 

CH2,2I-C1I  i»l  I 

PRU09303 

97* 

CK2, 31:0(1,31 

PRUO*  900 

98* 

Cl ( 3 t 1 ) :C1 ( 1 • 3 ) 

PRU09500 

99* 

C1(3,2I:C1<1,3> 

PRU09600 

SO* 

C 

CALCULATE  CSUBIJ»S  for  matrix 

PRU09700 

51* 

NUZX(2)=(EX(2I*NUXZ(2II/EZI2) 

PRU09833 

52* 

0ELTA  =  1.  .  NUX Y ( 2 1 -  Z.«NUXZ(2)*NUZX(2)  -  2. 

•N JXY ( 21 *NUXZ < 2 1 

PRU0*90Q 

53* 

.  *NUZX<2| 

PRU05000 

59* 

C2(l.l>=EX(2I*fl.  -  Nuxz<2»*MUZXf2n/DELTA 

PRU05100 

55* 

C2( l»2)S£X(2l*(NUXY(2i  «  NUX Z 1 2 ) *NuZ XI 2) l/OELTA 

PRU05200 

56* 

C2 ( I >3 1 =EX ( 2 I* ( I ,  *  NUXY(2 1 )«NUXZ (2 i /OELTA 

PRU05300 

57* 

C2(3,3>:EZ(2)*(1.  -  NUX Y 1 2  I **2 I /OELTA 

PRU05900 
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58* 

C2I2.1)=C2I1,21 

PRUOSSOO 

59* 

C2l2,2t=C2»l,l» 

PRU05600 

60* 

C2|2,3t=C2tl,3> 

PRU0S700 

6i* 

C2I  3,1>=C2I 1,31 

PRUOSSOO 

62* 

C2<3,2i:C2<l,3l 

PRUOj  900 

63* 

c 

CALCULATE  THERMAL  STRESS  COEFFS 

PRU07900 

66* 

BETAK  2 1  =BETA  x  HI 

PRU08000 

65* 

B£TA2|2|=SETA2«1I 

PRU08100 

66* 

00  30  1=1.3 

PRU08200 

67* 

gaHhahi»=0. 

PRUQb  300 

6S* 

00  30  JJ=1,3 

69* 

30 

6AMM*HI  I=6AHMA1II»  «  C  1 1 1  *  JJ»*BET  A 1 1  JJI 

70* 

00  60  1=1.3 

71* 

6AMMA2 (II -0* 

PRUqSTqo 

72* 

00  60  UJ?1,3 

73* 

60 

GAMhA2<II=GAHMA2IH  •  C2II.JJ)*3ETA2(JJ) 

76* 

c 

CALCULATE  ELASTIC  ANO  ThERnAL  STRESS  COEFFS  FOR  CONROSITE  AS 

PRU09200 

75* 

c 

FUNCTION  OF  FIBER  VOLUME  FRACTION.  Cll  *  C12.  C13,  C33.  GAMMA  1 , 

PRU09300 

76* 

c 

6AMMA3 

PRU09600 

77* 

Nl=VFt J) 

78* 

CALL  H| X 

PRUQ9800 

79* 

c 

CALCULATE  THERMAL  EXPANSION  COEFFS  FOR  COMPOSITE 

PRU09900 

80* 

BETAC(li:.S*l6AHMAE(l>*CEt3i3l  -  6ANNAE < 3 l*cC 1 1 . 3 t 1 /I .S«kE 

PRUlQOOO 

81* 

.  *CEI  3.3)  -  CE 1 1 • 3 ) **2 1 

PRU10100 

82* 

BETAEI3>:| .S*KE*GAMMAEI 3)  -  CE t 1 ,3» *6AMMAE | l » »/| .5*«i 

PRU10200 

83* 

.  *CEI  3.3)  -  CE<1.3>**2! 

PRU10300 

86* 

BETAS  1 21 =8ETAE < 1 ) 

PRU10600 

85* 

SAVEX:8ETAEt It 

86* 

SAVE2=BETAEI 31 

87* 

CINC7, jI=GAhhAEI 11 

88* 

CINIs. J»=S»MM4Et  31 

89* 

c 

THERM  CONO  FOR  COMPOSITE 

PRU10700 

90* 

N2CSN1«M2III  ♦  (1,  -  N j  !«KZ( 21 

pRUloBoo 

91* 

0=IKx<2>  *  XXI 11 1/IKX|1 1  «  KX(Z!  ♦  N1*IKXI2I  -  KXI1III 

PRU10900 

92* 

KXC=N1*mX(1)*I 1.  •  (1.  -Nl»*OI  ♦  tl.  -  Nl)*«X|2l*ll.  -  N1*QI 

PRUUODO 

93* 

c 

COMPOSITE  DENSITY 

PRU11100 

96* 

RH0C:N1*RH0I)|  •  (1.  -  N11*RH0I2I 

PRU11200 

95* 

c 

COMPUTE  Cll  ANO  C12  FOR  COMPOSITE 

PRU116O0 

96* 

CALL  MIXCC 

PRU11S00 

97* 

CALL  INVERT 

PRUH600 

98* 

CINtl.JI-CC(l.l) 

99* 

ClNI 2 , J 1 =CE ( 1.21 

300* 

C1N|3.J|=CEI1,3) 

101* 

CIN«6, jI=CEC 3,31 

102* 

c 

COMPOSITE  INPLANE  SHEAR  MODULUS 

“RU12200 

l(j3* 

C1N(5,J):6I2I*<2.*N1*g<1>  ♦  <1.  -  Nl»*«Glll  ♦  61 2 til / 1 2. *N1 *G|2 > 

106* 

.  •  11.  -  Ni l*IGl 1 1  ♦  GI2I1I 

105* 

ClNt6,J|=.S*(CEIl,ll  -  CE  <1,211 

106* 

c 

CALCULATE  COMPOSITE  ENGINEERING  MODULI  ANO  SAVE  F  DR  “TINTING 

PRU12600 

107* 

Ec*C=l./SI 1,11 

108* 

EC2C-1./SI 3,31 

109* 

NUCXTCr  -  S(l,2l*£CXC 

no* 

NUCX2C:  -  SI  1 , 31 *£CZC 

IU* 

NUCZXC=NUCX2C*CCXC/EC2C 

112* 

CALL  STRESSISIG21,SIGZ2,SI6R,SIGHll,sIGH21,SIGH2Z.SAVEX.SAV£ZI 

113* 

WRITE  1 6. 1o3I 

U6* 

WRITCI6.100)  J 

115* 

too 

FORMAT!/, *  PROPERTIES  OF  PLY',12,'  AS  CALCULATED  VIM  “TUFC'I 

U6* 

WRITE  16, 1011  £CXC,£CZC,nUCXVC,NUCXZC 

117* 

tot 

FORMAT  1  1H  ,SX,’ET:*,E11.S,5X,*EA  =  *,E11.S,SX,*NUTT=* ,rB.S.SX,«NUTA 

z 

118* 

•*.F8.S| 

H9* 

WRITE  16, 1021  CINIS, J| ,CIN(6,JI,BETAEI1I,BETAE!3I 

120* 

132 

FORMAT  Ih  ,SX,*GAT=*,Ell.S.Sxf*STT:*,Ell.S,SX,*ALFT:«,Ell.S,SX,*AL 

I2i  • 

.FA:', Ell. 51 

122* 

103 

FORMAT!  //I 

123* 

RETURN 

126* 

ENO 

PRU2Z603 

c»o  0 f 

COMPILATION:  NO  DIAGNOSTICS. 
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1* 

SUBROUTINE  ROOTlN.XOUTI 

R0000100 

2* 

C  SOLVES  N  LINEAR  SIMULTANEOUS  EQUATIONS.  N.LE.6  UNLESS 

R0000200 

3* 

C  DIMENSIONS  ARE  CHAN6E0 

R0000300 

«• 

IMPLICIT  OOUsLC  PRECISION! A-H.O-Z) 

R000090Q 

S* 

REAL  X0UTI6),  XKINI 6, 7> 

ROOOOSOO 

6* 

DIMENSION  X ( 6  }  .  XK | b.7) 

R0000600 

7* 

commOn/coef/xkin 

R0000700 

S* 

NP  I - N  •  X 

Roooosoa 

9* 

NM I 5N  -  1 

R0000900 

0« 

NI  =  0 

ROOOXOOO 

X* 

00  ID  1  =  1, N 

ROOOXXOO 

2* 

DO  ID  JsX.NPI 

R000X20Q 

3* 

XO  XKlI,JI=XKINII,J) 

RO 00 1300 

*« 

20  N1=NI  ♦  X 

ROOD! 900 

5* 

IMNI.GT.NMXI  GO  TO  80 

RO  001 500 

6* 

NIPl=iI  ♦  I 

R0001600 

1* 

A;XK|N1,NII 

R0O0X700 

8* 

IM-NI 

R000180Q 

9* 

Do  30  I=NIP1«N 

ROOOX  900 

0* 

IFIOABS 1 XKI I ,NI 1 1 .LE.OABSI Al 1  GO  TO  30 

R0002000 

1* 

A=XKII,NI) 

R0002100 

2* 

IM=I 

90002200 

3« 

30  CONTINUE 

R0002300 

9* 

lFHM.LE.NIt  GO  TO  50 

R0002900 

s* 

00  90  J=NI,NP1 

R0002500 

6* 

AXXvHh.JI 

R0002600 

7* 

XNI1M,jI=XN(NI.JI 

R0002703 

8* 

«0  XKINI,Jt:A 

R0002800 

9* 

SO  CONTINUE 

RO  002900 

0* 

a:xk(N1.nii 

R0003003 

X* 

00  60  J=NI »NP 1 

R0003X00 

2« 

60  XK(Nl,J)=XK<NI,JI/A 

R0003200 

3* 

00  70  IJTNIPX.N 

R0003300 

9* 

00  70  JSNIPX.NPI 

R000390Q 

5* 

70  XK<I.Ui:XNf I.NltaXHINI.JI  -  XKII.Jt 

R0003503 

6* 

GO  TO  20 

R0003600 

7* 

80  CONTINUE 

R0003700 

a* 

XI Nt=XKI N.NP1 I/XKIN.NI 

R0003800 

9* 

NX  =  N 

ROOq39qq 

0* 

90  NX -NX  -  1 

R0009003 

X* 

NXP1-NX  ♦  I 

R0009X00 

2* 

IFtNX.LT. XI  60  TO  1X0 

R000920Q 

3* 

X 1  NX )=XKtNX, NP1 1 

R0009300 

9* 

00  100  MTNXPl.N 

R0009900 

5* 

XOq  X|NX):XINXI  -  XK (NX ,K )*Xf N  ) 

ROOQ9S33 

8* 

X(NXI=XINXI/XKINX,NXI 

90009  600 

7* 

GO  TO  90 

R0009700 

8* 

no  continue 

R0009  BOO 

9* 

00  X20  1=1, N 

RO  009  9qQ 

0* 

120  XOUT ( I 1 :xf I 1 

ROOOSOQ3 

X* 

RETURN 

R0005100 

2* 

ENO 

R0005200 

N3  OF  COMP  I L  A  T I  ON  t  NO  DIAGNOSTICS* 
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1* 

2* 

3* 

9* 

5* 

6* 

7* 

8* 

9* 

10* 

11* 

12* 

13* 

18* 

IS* 

18* 

17* 

18* 

19* 

20* 

21* 

22* 

23* 

28* 

25* 

26* 

27* 

28* 

29* 

30* 

31* 

32* 

33* 

38* 

35* 

3b* 

37* 

38* 

39* 

*0* 

*1* 

82* 

83* 

99* 

95* 

96* 

87* 

98* 

89* 

SO* 

51* 

52* 

53* 

59* 

55* 

56* 

57* 

58* 

59* 

60* 

61* 

62* 

63* 

68* 

65* 

66* 

67* 

69- 

69* 

70* 

71* 

72* 

73* 

79* 

75* 

76* 

77* 


SUBROUTINE  SCTUPINI 

C  SETS  UP  MATRIX  OF  COEFFS  USED  B»  SUBROUTINE  NIXCC  To  03TAIN 

C  Cll  ANO  C12 

REAL  Nl,  M I  6, 7 1 

dimension  etai2i,  C66i2i 

COMMON/EL  AST /Nl ,  c  IH21,  C12I2I 

COMMON/COEF/R 

IFIN.EO.zl  60  TO  10 

R1=S0R7INH 

R12:R1**2 

R18=R12**2 

R16=R12**3 

etai  n=ci2(  ii/cimi 

ETAI2I=C12<2I/C11I2I 

C66|ll=.5*lCllf II  -  C 12  ( 1 1 1 

C66I21  =  .5*IC1K2I  -  C12I211 

NI1,11=R19 

Ml I,2»=RI6 

K11,3I=R12 

Ml  1,81:1. 

Mli,5»=-Bl9 
Ml 1,61=-R16 
Ml l, 71=0. 

MI2,1I=-£TAI jl*ETAI2)*R19 
MI2,2»:-.5*ETAI1 1*13.  *  ETAI2I  )*R16 
Kl2,31=-.5*£TAlll*CTAI2)*ll.  -  ETA|2»I«R12 
MI2,*l  =  ET»m*CtAI2» 

KI2,5»=£TAI1»*ETAI2J*R18 
K|2,6I=.5*ETA|21*I3.  •  ETAI|I)*R1b 
9(2,71=0. 

MI3,1I=R19*IC11(21  -  Cl  2  I  2  1  1 
Ml  3,21  =  0. 

Ml  3,3>  =  -Ri2*ICUI2>  -  ETAI2I*C12I2»» 

MI3,8»=-3.*|C11|21  *  Cl 2 1 2 1 1 
K(3,51  =  -Rl9«tCllUI  -  C12I1I1 
Ml  3,61=0. 

M I  3, 71  =  0. 

M,8,1|=-2.*ETAIII*£TAI2»*Rj%*C6bI2» 

KI9,2l=-C66(2>*R16*ETAt ll*3.*ll.  ♦  ETAI2II 

MI8, 31=-C66I2»*ETAI1I*CTAI2»*I1.  ♦  EVAI2II*R12 

Ml9,*l=-6.*C66I21*ET*III*ETAI21 

M|8,5»=2.*C66I11*CTAUI*CTAI2»*RI8 

919, 61  =  3. *C66ll»*ETA|2l*ll.  «  ET  A  1 1  I  I  *R  16 

KIM, 71=0. 

MI5,  11  =  C11«2»  -  C12 1 2 ) 

MIS, 21=0, 

MI5,3|=~C11I21  ♦  ETAI21*C12l2l 
K(5,81  =  -3.*(CHI21  -  C12I211  ’ 

MIS,51=0. 

MIS, 61=0. 

K|5,7I=1. 

KI6,ll=-2.*CTAI2)*C66|2l 
MI6,21=-3.*C66I2»*I 1.  ♦  Et*I 2  1 1 
M I  6, 31  =  “C66I 2I*ETaI2I*I 1.  ♦  ETAI211 
M|6, 91= •6,*ETAIjI*C 6^121 
Ml  6, 51  =  0. 

K 16,61=0. 

MI6,71=-CTAI21 
RETURN 
10  CONTINUE 

ETA|21=C12I21/C11I2I 
MIS, 11=1. 

K (5 ,2  >  =1 . 

M I S , 31  =  1 . 

MIS, 91=1. 

K I  5 , 51  =  0 • 

KlS,6l=0. 

kIS,7I=1. 

K I  6 , 1 1 =  -  CTAI21 

N 1 6,21 =  -  .5*13.  ♦  ETAI2II 

M 1 6 , 3 1 =  -  .5*11.  -  ETAI 2 1  I *E  T  n I  2  ) 

K I  6, • I =ETA 1 2 1 

K 1 6 , SI =0* 

K 16,6 1 =0. 

K 1 6 , 71 =  -  ETA  121 

RETURN 

ENO 


SETQ0200 
S£ TOO 300 
SETQ0900 
SETQOSOO 
SE7Q0633 
SET00700 

SET00800 

SET00900 

SE  TQ) 000 

SETOUOO 

SET01203 

SETQ1300 

SET0190Q 

SET01500 

SETQ16Q0 

SE 701700 

SETQ1803, 

SE  TQ1900 

SET02000 

SE  T02100 

SE  T022DD 

SET02300 

SE  T02900 

SE  702530 

SE 702600 

SET02700 

SETQ2800 

SET02900 

SETQ3000 

SET03100 

SET03203 

S£  T03300 

SETQ3900 

SET03500 

SE T03600 

SE  TQ37qq 

SE  TQ383Q 

SE  Tg3900 

SE 709000  I 

sEtomioo 

SE 709200 
SCTQ9300 
SE 708830 
SE 708500 
SE  T0960Q 
SET09700 
SE  T09 8qo 
SET09900 
SETqSOOO 
SE  705100 
SETQ5200 
SE  705  300 
SET05900 
SC  705500 
SE  705600 
SC  TQs  700 


SET0560Q 

SCT0590C 
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2 

2 


C 


SUBROUTINE  STRESS ISISZ1 ,SIGZ2,SI6R,SIGHll,SISM2l,Sl3422,8CTAX, 
•  BETA Z ) 

CALCULATES  INTERIOR  FIBER  ANO  MATRIX  STRESSES 
REAL  N1 

COMNON/MIXl/Clt 3.31 .  C2  (  3 . 3  I  .  6AMMA I  131.  GAMMA  2131.  3ET41I3I. 

.  BE  T  A  2 1  3  I 

COMMON /EL AST /N 1 .  C1H2I,  C12I2I 

A2-IC2I  1.3  )*BETA2  ♦  IC2II.il  ♦  C2U.2I  IRBETAX  -  6AMm42I1I)/ 

.  INl«2.«C2(l,ln 
B2=BETAX  -  A2*NI 
B1-A2  •  B2 

SIG2l=2.*CH  1,31*81  ♦  C1I3,3>«BETA2  -  6AnmAH3> 

SI622=2.*C2l 1,3>*B2  ♦  C 2 1  3 , j 1 «BE T A2  -  6AMMA2I3I 
SI6r:IC1I1.1I  «  CHl,2n*Bl  ♦  C1U.3IRBETA2  .  6AMMA 1111 
SlGMll=SlGR 

SI6H21=C2< 1.2I»<-a2  ♦  821  ♦  C2d,ll*<A2  ♦  321  ♦  C2t  1,3I*8ETA2 
.  -  G A  HMA  2  111 

S16H22=C2(  1,2)*(~A2«N1  ♦  821  ♦  C2 ( 1 , 1 > • ( A 2*M  1  ♦  821  »  C2I1.3I 
.  «BET AZ  -  GAMMA 2  I  II 
RETURN 
ENO 


STROOlOa 
STR00200 
"  STR00300 
strooaoq 
stroosoo 
STR00600 
STR00700 
STROOSOO 
STR00900 
STROIOOO 
STROUqo 
STR01200 
STRq13oO 
STROIAOO 
STR01SOO 
S  TR01 633 
STROl TOO 
STRQ1SDQ 
STROl 9OO 
STRQ2000 
STR02100 


I  OF  COMPILATION: 


NO  DIAGNOSTICS. 
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1* 

2* 

3* 

9* 

5* 

l» 

7* 

I* 

8* 

10* 

11* 

12* 

13* 

1** 

15* 

16* 

17* 

18* 

19* 

23* 
21* 
22* 
23* 
29* 
25* 
2b* 
27* 
28* 
29* 
30- 
SI* 
32* 
33* 
39* 
35* 
3b* 
37* 
38* 
39* 
93* 
91* 
92* 
93* 
99* 
95* 
9b* 
97* 
98* 
99* 
50* 
SI* 
52* 
53* 
59* 
55* 
Sb* 
5  7* 
58* 
59* 
b0* 
bl* 
82* 
63* 
69* 
65* 


SUBROUTINE  ThERnlIAlF  •SISOOTI 
0 1  ME  NS I  ON  Xlfc),  ALF 19  1,  SI6oOT(9,3l 
REAL  Nl,  N2«  N3,  N 
C ONHON/ ONE /N 1 ,  N2,  N3 

COhhoN/TU0/CH9,*I,  C2 I  9 ,9  )  ,  0319,9) 

CONHON/gOEF /* I b, 7 ) 

C0MN0N/THREE/GANM*1I*I,  5ANN«2(9),  6AHMA3I9I 
*11,11=01(1,1) 

8*1,21=0. 

*11,31=0. 

*<1,91=0111,21 

*(i,si=cm,3i 

*11, 61=2. *0(1,9) 

Nll,7l=6AHNAxlll 

*12,11=0. 

*12,2 1=C2I 1,11 
N I  2, 31=0. 

*12,91=021 1,21 
*12,51=021 1,31 
*I2,6I:2.«C2|1(9I 
*12,7 1 =GAHMA  2111 
*13,11=0. 

*13,21=0. 

*13,31=0311,11 
*13,91=031 1,21 
*13,51=031 1,31 
*1 3,6l=2.*C3t 1,91 

*  1 3, 71  =  6AHNA31 1 1 
*19, 11:91*01 (1,21 
*I9,21=n2*C2I 1,21 
*19, 3t=N 3*0311,21 

*I9,9I=N1*C1I2,21  ♦  N2* 0212, 21  *  n3*C3I2,21 
KI9,5l=Nl*Cil2,3l  ♦  N2*C2 12,31  ♦  N3*c312,3l 

*  19, bl:2.*(Nl *0112,91  ♦  N2*C2 1  2 ,9 1  ♦  93*03(2.9)1 
*19,71 =Nl*GAHMA 112)  ♦  92*6 ANN A  2 ( 2  )  .  N3*GANNA  3  1 2  I 
*15, 11 =N 1*0 111, 3) 

*15,2) =N2*C2I 1,3) 

*  ,5, 31=93*03 11,3) 

H15,9)=N1*C1|2,3)  ♦  92*0212,3)  »  93*03(2.3) 

*IS,S)=n1*C1I3,3I  •  N2*C2 13,3)  «  N3*C3l3,3) 

*  I  5, 61  =  2. *191*01 1  3,9)  ♦  n2*C2I 3,9 1  *  N3*C3I3,91) 
KI5,7)=9i*6ANNA1I3)  «  N2*GANN«2I3>  *  N3*GANNA3I3) 
*lb,l)=Nl*Clll,9l 

Klb,2)=N 2*0211,91 
*lb, 31=93*0311, 91 

*16,91=91*0112,9)  «  92*0212,9)  ♦  *3*0312.91 

Klb, 51=91*0113, 9)  ♦  N2*C213,9I  •  93*0313,9) 

*lb,bl:2. •(91*0119,9)  «  92*0219,9)  .  93*0319.9)1 
Mlb,7)  =  91«GANMAll9)  •  N2«6ANNA2(9)  ♦  93*6 ANNA  3(9) 
CALL  ROOTlb.X) 

EPSX  =  9 1«X  111  «  9  2*X I  2 1  *  9  3  * X  I  3  I 
*LF 1 1 1 =ERSX 
ALF (2|=X|9) 

ALFI 3)=XI5) 

ALf l9)=Xlb) 

SIGOOTI 1,11=011 1,1I*XI1  1  •  Cl  I  1 , 2  1  *X I  9  I 
.  *  2.«Clll,9)«Xtbl  -  GANNA  111) 

SIGDOTI2,l)=Cltl,2l*X(l)  *  C1I2,2)«X(9) 

.  ♦  2.*C112,9)*Xib)  -  GANNA1I2) 

SIGOOTI3, 11=0111, 3)«XI1  )  •  Cl  I  2, 3  )  *X 1 9  ) 

.  ♦  2 ,*C 1 1  3 ,9 1 • X 1 6  »  -  GANNj 1 1 3  I 

SIGOOTI9, 11=0111, 9>*XU  )  ♦  C)I2,9)*X(9) 

.  •  2,«Cll9,9)«Xtb)  -  GANNA1 (9  | 

SlGOOTIl, 21=0211, 1)*XI2)  *  C2tl,2)*XI9l 
.  ♦  2.*C2*1,9)*X(6)  -  G»NN»2 1  1 ) 


♦  0111,31**15) 

♦  Cl  12,3 ) •* I S  ) 

♦  0113, 31 *X 15) 
«  CH  3,9)  *XI  5) 

♦  C2(1,3)*XI5) 
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66*  Sl600TI2,2«=C2«i.2)*XI2»  ♦  C2«2,2>*Xt»»l  ♦  cZ«2»3»**«5 

67*  •  *  2 • *C2 1  2  »* )«X (6 l  -  GAMMA2I2I 

66*  SIGOOTI  3.2>=C2<  1  ,3)*X  C2  I  ♦  C2  C  2 , 3  >  *X  I  •*  |  ♦  C2t3,3»*XIS 

69*  .  ♦  2.*C2< 3.*>*X(6>  -  GAMMA2«3» 

70*  SIG00TM,2ISC2»1.6»*XI21  ♦  C2i  2  »*  J  *X  I  <•  I  ♦  C2<3,«>«XIS 

71*  .  ♦  2 . «C2 1* ,6 )*X 16  I  -  GAMMA2I9I 

72*  SIG00r(l«3)=C3ri(l)«XI3)  *  C3(t,2>*X|»)  *  C3U,3»*XI5 

73*  .  «  Z.*C3(l'«>*Xtb>  -  GAHMA3(1> 

79*  SIGDOTT 2»  3  I - C3* l »  2>  *X  <  3  I  ♦  C3« 2. 2 > *X 1 9 >  *  C3»2,3»«x<5 

75*  .  •  2 . *C3  I  2  ,  H  )  *X ( b  I  -  GAMMA  312) 

76*  SlGDOT(3,3):C3ll,3l*XC3l  ♦  C3( 2. 3 1 *X t « I  ♦  C3»3.3>*XI5 

77*  .  *  2.*C3« 3.4»**»6>  -  GANMA3 4  31 

78*  SIGDOT<9«3>-C3<ltHl*X(3l  ♦  C  31  2  » *  I  *X  I  9  I  •  C313.9»*X«5 

79*  .  ♦  2.*C3l9f9»*X«6 J  -  GAMMA 3 18  1 

80*  RETURN 

81*  ENd 


END  OF  COMPILATION*  no  diagnostics. 
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